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1.  Introduction 


FASCODE  for  the  Environment  (FASE)  is  the  latest  in  a  long  line  of  atmospheric 
transmittance  and  radiance  models  developed  by  the  Optical  Physics  Division  of  the  Geophysics 
Directorate,  USAF  Phillips  Lab  (Anderson  et  al.,  1995).  Of  interest  to  this  project  is  the  high 
spectral  resolution  Fast  Atmospheric  Signature  Code,  “FASCODE”  (Smith  et  al.,  1978).  The 
most  recent  version,  FASCOD3P,  was  released  circa  1992.  Line-by-line  models  such  as 
FASCODE  calculate  atmospheric  molecular  transmittance  from  first  principles,  the  core  of 
which  is  a  fast  algorithm  to  calculate  the  Voigt  line  shape.  The  Voigt  algorithm,  coupled  with 
the  best  available  physics  and  latest  atmospheric  models  and  data,  makes  FASCODE  an  accurate 
and  efficient  model  for  atmospheric  transmittance  and  radiance.  Further,  it  is  accepted  as  a 
reference  standard  model  for  the  scientific  community. 

After  the  FASCOD3P  release,  work  on  the  model  continued  at  AER  in  two  different 
directions.  Internally,  AER  modified  FASCOD3P  for  application  to  inversion  problems  and 
sensor  studies  (this  code  is  called  XFWD  and  is  described  in  Miller  et  al.,  1995).  The 
Department  of  Energy  (DOE)  Atmospheric  Radiation  Measurement  Program  (DOE-ARM)  also 
sponsored  AER  for  the  development  of  the  Line-By-Line  Radiative  Transfer  Model  (LBLRTM), 
with  the  main  goal  of  further  improving  the  physics  though  the  analysis  of  atmospheric 
measurements  made  at  various  ARM  experimental  sites  (Clough  et  al.,  1993).  LBLRTM  is 
derived  from  FASCOD3P  and  includes  a  number  of  significant  advancements  and  changes. 

Both  XFWD  and  LBLRTM  are  significantly  different  in  capabilities  from  FASCOD3P. 
In  particular,  certain  capabilities  of  FASCOD3P  (multiple  scattering,  non-LTE  emission,  and 
aerosol  attenuation)  were  removed  for  these  specialized  applications.  The  FASE  model  was 

developed  to  continue  the  FASCODE-series  of  models  while  meeting  the  following  criteria: 

-  incorporate  the  best  features  and  advancements  from  LBLRTM  and  other  models, 

-  preserve  all  the  capabilities  of  FASCOD3P, 

-  improve  the  user  interface,  and 

-  add  other  features  and  improvements  as  resources  permit. 

This  report  describes  the  work  performed  in  developing  FASE  and  addresses  upgrades  to 
existing  FASCOD3P  features,  the  addition  of  new  features,  numerical,  algorithm  and  coding 
comparisons  among  the  various  models,  and  suggestions  for  future  code  development. 
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2.  Model  Comparisons 


In  choosing  the  starting  point  from  which  to  develop  FASE  we  compared  three  line-by-line 
radiative  transfer  codes:  FASCOD3P  (Phillips  Laboratory),  a  modified  version  of  FASCOD3  (an 
internal  product  of  AER  not  intended  for  public  release),  and  the  "global  release",  September 
1994,  version  of  LBLRTM  (AER  /  ARM).  It  should  be  noted  that  all  three  models  are  similar  in 
that  they  were  derived  from  earlier  versions  of  FASCODE;  the  modified  version  of  FASCOD3 
can  be  considered  an  intermediate  step  between  FASCOD3  and  LBLRTM.  The  object  of  this 
internal  study  was  to  determine  which  of  the  three  models  (FASCOD3P,  the  modified  version  of 
FASCOD3,  or  LBLRTM)  would  allow  us  to  spend  a  minimal  amount  of  time  re-coding  features 
from  other  models  so  that  our  time  could  be  spent  updating,  and  adding  to,  the  model  physics. 
This  model  comparison  was  accomplished  through  a  side-by-side  examination  of  the  three 
models.  In  addition  we  also  investigated  the  cause  of  observed  differences  in  the  radiances 
calculated  with  F ASCOD3P  and  with  LBLRTM,  including  both  timing  and  numerical 
comparisons.  The  following  sections  report  the  results  of  these  comparisons. 

Section  2. 1  describes  the  differences  in  algorithm  features  and,  where  applicable,  their  impact 
on  the  accuracy  of  the  calculation;  Section  2.2  describes  numerical  differences  between  the  codes 
for  a  set  of  test  cases. 

2.1.  Features  Comparison 

Table  1  is  a  partial  list  of  the  distinguishing  characteristics  of  the  three  models.  It  provides  a 
comparison  of  FASCOD3P,  the  AER-modified  version  of  FASCOD3,  and  a  "global  release" 
version  of  LBLRTM.  It  should  be  noted  that  while  this  is  only  a  partial  list  of  the  differences,  it 
represents  features  important  to  the  development  of  FASE. 
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Table  1  Model  Comparisons 


Model 

Features 

FASCODE3P 

Modified 

FASCOD3 

LBLRTM 

Voigt  Function 

Old 

Old 

New 

Vectorized 

No 

No 

Yes 

Parameterized 

None 

Some 

Extensive 

Line  Rejection  Flag 

No 

No 

Yes 

Multiple  Scattering 

Internal 

None 

Output  in  correct 

LASER  Option 

Yes 

No 

format  for  external 
codes  (e.g.,  CHARTS 
or  DISORT) 

No 

Non-LTE  Option 

Yes 

No 

No 

Scanning  Function 

SCANFN 

SCANFN 

FFTSCAN  and 

Plotting  Function 

FASPLT 

FASPLT 

SCANFN 

LBLPLT  with  NCAR 

Aerosol/Cloud 

Yes 

No 

Graphics  Option 

Yes 

Attenuation 

Spectral  Bandwidth 

520  cm'1 

520  cm1 

2020  cm'1 

Limit 

I/O  Routines 

Old 

New 

New 

(BUFIN/BUFOUT) 

Function 

Geometry  Package 

FSCATM 

FSCATM 

LB L ATM 

2.1.1.  Voigt  Function  Algorithm 

Atmospheric  infrared  absorption  lines  are  described  by  the  Voigt  line  shape,  which  is  the 
convolution  of  the  Lorentz  and  Doppler  line  shapes.  The  Voigt  function  has  no  analytic 
expression  and  considerable  effort  by  many  authors  has  been  devoted  to  finding  efficient  and 
accurate  approximations.  At  the  core  of  the  FASCODE  model  is  a  very  fast  Voigt  algorithm  that 
calculates  the  Voigt  function  via  a  lineshape  decomposition  (Clough  and  Kneizys,  1979). 
However,  recent  advances  in  the  accuracy  of  field  measurements  and  the  corresponding  more 
stringent  requirements  for  simulations  require  even  greater  accuracy  for  the  Voigt  function. 
LBLRTM  has  modified  the  Voigt  function  from  FASCODE  to  increase  the  accuracy.  This  work 
is  one  of  the  driving  factors  behind  FASE  as  it  could  significantly  impact  the  calculations  under 
certain  scenarios,  particularly  atmospheric  parameter  retrieval  problems. 
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The  FASCODE  lineshape  algorithm  decomposes  the  Voigt  lineshape  function  into  a 
summation  of  functions  representing  the  Lorentz  and  Doppler  components  as  well  as  an  error 
term  to  account  for  small  differences  from  these  functional  approximations: 


V£,z)=  C,  (E)F|  (z) + C2  (5)F2  (z) + Cj  (z) 

+-^[F,(z)  +Fs(z)]+Cd(^d(z) 
■«v,(5)V,(z) 

where 


Eq.  1 


a 

ac  +  «b 


Eq.  2 


and 


V-U 
z  - - 1 


a 


v 


In  these  equations: 


Eq.  3 


v  =  frequency  in  wavenumbers 

Vj  =  frequency  of  the  i’th  spectral  transition 
av,ac,andaD  =  the  Voigt,  Lorentz  and  Doppler  halfwidths 
Cj,  i  =  l,2,3,D,e  =  coefficients 
F.,  i  =  1,2,..5,D  =  shape  functions 
Ve  =  error  shape  function 


The  Lorentz  component  is  described  by  the  first  four  functions  in  the  Doppler  with  a  single 
function  (subscript  D),  and  the  error  with  the  remaining  term.  The  coefficients,  C(£),  are 
tabulated  in  data  statements  within  the  code,  and  the  shape  functions,  F(z),  are  computed  within 
the  code. 


In  both  FASCODE  and  FASE  there  are  102  coefficients  stored  in  data  statements.  However, 
to  determine  the  coefficient  corresponding  to  the  value  of  FASE  linearly  interpolates  between 
coefficients  while  FASCODE  merely  chooses  the  coefficient  closest  to  the  calculated  value  of 
The  shape  functions  are  also  computed  in  a  slightly  different  manner:  FASE  computes  F(z)  on  a 
grid  of  2001  points  while  FASCODE  computes  only  201  points.  The  combination  of  these  two 
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changes  in  FASE  decreased  the  errors  associated  with  the  lineshape  computation  as  compared  to 
measurements  (Clough  and  Brown,  1994),  but  increased  the  overall  computation  time.  In  FASE 
this  timing  penalty  was  offset  through  the  optimization  of  other  parts  of  the  algorithm  coding, 
namely  the  BUFIN/BUFOUT  routines.  This  study  evaluates  the  differences  between  the  original 
FASCODE  method  and  the  FASE  improvements  in  order  to: 

•  determine  whether  or  not  the  changes  are  justified  in  terms  of  accuracy  /  timing 
considerations,  and 

•  decide  whether  there  is  a  less  computationally  expensive  way  to  achieve  the  numerical 
accuracy  of  the  increased  points  and  interpolation. 

One  should  note  that  the  number  of  points  and  interpolation  differences  appear  in  the 
convolution  routine  for  the  first  three  functions,  while  the  so-called  fourth-function  is  different 
because  it  only  interpolates  the  halfwidth  ratio  (the  same  number  of  points  are  used  to  define  the 
lineshape  function).  Because  the  interpolation  time  for  the  fourth-function  is  very  small 
compared  to  the  other  functions,  comparison  of  fourth-function  timing  gives  a  measure  of  the 
numerical  stability  of  the  timing  test. 

2.1.1 .1.  Voigt  Function  Timing  Tests 

The  comparison  of  model  timing  can  be  difficult  since  the  time  will  vary  depending  on 
machine  load  and  network  traffic.  To  minimize  this  effect,  30  runs  of  the  models  were 
compared.  Further,  the  runs  were  done  after  midnight  on  a  (nearly)  dedicated  machine  to  reduce 
the  competition  for  CPU  and  disk  access  time.  A  total  of  three  different  model  configurations 
were  used  to  ascertain  timing  differences:  FASE,  FASCODE,  and  FASE  with  only  201  points 
for  the  shape  functions,  hereafter  denoted  as  FASE2001,  FASCD3P,  and  FASE201.  Because  of 
coding  differences  such  as  BUFIN/BUFOUT  (the  file  I/O  routines),  the  layer  convolution  times, 
rather  than  the  total  run  time,  were  the  basis  for  the  comparison  (a  multiple-layer  profile  was 
used  in  order  to  increase  the  number  of  points  for  comparison).  Comparison  of  FASE2001  with 
FASCD3P  will  show  the  overall  change  in  CPU  time,  comparison  of  FASE2001  with  FASE201 
will  illustrate  the  additional  time  due  to  the  increased  number  of  points,  and  comparison  of 
FASCD3P  with  FASE201  will  show  the  effect  of  the  coefficient  interpolation.  To  minimize 
differences  in  the  overall  timing  of  the  codes,  only  the  optical  depth  was  computed. 

Three  types  of  tests  were  used  to  compare  the  convolution  times: 

1.  20  layer  atmosphere  from  0  to  25  km  (looking  up), 

2.  22  layer  path  from  100  km,  through  a  tangent  at  25  km,  to  space,  and 

3.  22  layer  path  from  50  km,  through  a  tangent  at  5  km,  and  back  up  to  50  km. 
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All  tests  were  over  the  spectral  interval  from  2050  -  2075  cm'1;  subtle  differences  between 
FASE  and  FASCODE  were  eliminated  by  explicitly  specifying  the  atmospheric  layering  and  by 
setting  ALFALO  to  0.04  cm  *.  The  tangent  path  of  Case  2  was  tested  twice,  once  after  midnight 
on  a  Sunday  morning,  to  judge  the  effect  of  computer  load  on  the  timing.  The  results  of  these 
tests  are  given  in  Table  2  through  Table  4,  summarizing  the  layer  convolution  times,  total  model 
run  times,  and  total  time  savings. 

The  percent  difference  in  layer  convolution  times  between  the  model  configurations  were 
computed  for  each  layer  and  averaged  over  all  of  the  runs.  "Mean"  refers  to  the  mean  ratio  over 
all  the  layers  and  all  of  the  test  runs,  "RMS"  is  the  root-mean-square  variation  from  this  mean, 
and  "Savings"  is  the  percent  time  saved,  equal  to  (1-Mean).  The  ratio  of  FASCD  to  FASE2001 
shows  the  effect  of  increasing  the  number  of  points  used  to  describe  the  line  shape  and 
interpolating  between  the  coefficients,  FASCD  /  FASE201  shows  the  effect  of  the  interpolation, 
and  FASE201  /  FASE2001  shows  the  effect  of  the  increased  number  of  points.  The  "savings" 
given  in  the  first  row  of  each  case  should  be  the  sum  of  the  next  two  rows,  but  this  is  not  exactly 
true  because  of  differences  between  other  portions  of  the  FASE  and  FASCODE  algorithms.  The 
ratio  of  FASE  runs  for  LBLF4  gives  an  indication  of  the  stability  of  the  test  since  this  subroutine 
is  the  same  for  both  of  the  FASE  models. 

For  the  up-looking  case  (Case  1)  and  for  the  lower-altitude  limb  case  (Case  3),  eliminating 
the  interpolation  saves  more  time  than  decreasing  the  number  of  points.  The  opposite  is  true  for 
the  high-altitude  limb  tests  (Case  2).  This  illustrates  that  determining  the  effect  of  changes  to  the 
model  involves  a  complex  relation  between  the  algorithm  itself  and  the  characteristics  of  the 
atmospheric  path.  One  should  also  note  that  in  the  high-altitude  tangent  cases  the  time  saved  by 
model  changes  is  much  less  than  for  the  other  test  cases.  This  could  be  due  to  the  higher  spectral 
dv  required  for  this  case.  In  any  event,  this  illustrates  the  difficulty  of  evaluating,  for  all 
scenarios,  the  timing  impact  of  model  changes. 
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Table  2  Percent  Difference  in  Layer  Convolution  Time  Between  Model  Configurations 


Case 

Ratio 

Mean 

HIRAC 

Savings 

RMS 

Mean 

LBLF4 

Savings 

RMS 

1 

FASCD/FASE2001 

81.69% 

18.31% 

5.37% 

105.71% 

-5.71% 

17.74% 

FAS  CD/F  ASE20 1 

86.35% 

13.65% 

5.67% 

105.74% 

-5.74% 

17.68% 

FASE  201/2001 

94.61% 

5.39% 

1.45% 

99.98% 

0.02% 

1.81% 

2A 

FASCD/FASE200 1 

94.36% 

5.64% 

98.44% 

1.56% 

3.87% 

FASCD/FASE201 

99.51% 

0.49% 

1.63% 

98.23% 

1.77% 

3.50% 

FASE  201/2001 

94.83% 

5.17% 

1.66% 

100.21% 

-0.21% 

1.31% 

2B 

FASCD/FASE2001 

94.15% 

5.85% 

98.29% 

1.71% 

3.35% 

FAS  CD/F  ASE201 

99.40% 

0.60% 

0.86% 

98.20% 

1.80% 

3.32% 

FASE  201/2001 

94.72% 

5.28% 

0.30% 

100.09% 

-0.09% 

0.09% 

3 

FAS  CD/FASE200 1 

87.75% 

12.25% 

5.86% 

102.25% 

-2.25% 

8.97% 

FASCD/FASE20 1 

92.98% 

7.02% 

5.70% 

102.62% 

-2.62% 

8.87% 

FASE  201/2001 

94.37% 

5.63% 

2.41% 

99.65% 

0.35% 

2.38% 

Table  3  Total  Model  Run  Time  Comparisons 


Table  4  Evaluation  of  Time  Saved 


Case  1 

Case  2A 

Case2B 

Case  3 

HIRAC 

Savings  without  interpolation 

4.51% 

0.03% 

0.04% 

1.10% 

Savings  without  increased  points 

1.78% 

0.30% 

0.31% 

0.88% 

Total  Savings 

6.04% 

0.33% 

0.35% 

1.92% 

Time  Saved 

6.72  sec 

0.98  sec 

1.01  sec 

3.83  sec 

TOTAL 

Savings  without  interpolation 

3.55% 

1.18% 

1.24% 

-0.21% 

Savings  without  increased  points 

6.05% 

0.16% 

0.25% 

1.06% 

Total  Savings 

5.09% 

1.34% 

1.50% 

0.80% 

Time  Saved 

5.67  sec 

3.98  sec 

4.32  sec 

1.60  sec 

A  comparison  of  total  model  run  times,  expressed  in  seconds,  is  given  in  Table  3.  The  mean 
value  is  over  the  30  runs  tested  for  each  case.  Case  2B  was  run  after  midnight  on  a  Sunday 
morning  and  the  RMS  variations  in  the  total  times  are  much  less  than  for  the  other  cases.  Note 
that  the  total  run  times  are  not  indicative  of  the  timing  differences  between  the  models  because  of 
other  differences  in  the  algorithm  coding.  Case  2  was  significantly  longer  than  the  other  cases 
because  the  higher  altitudes  require  calculations  on  a  smaller  spectral  grid. 

Given  the  relative  changes  in  the  individual  convolution  routines  and  the  total  time  to  run 
each  model,  it  is  straightforward  to  calculate  the  amount  of  time  that  will  be  saved  by  making  the 
changes.  However,  the  models  have  slightly  different  coding  for  many  of  the  routines,  due  in 
part  to  the  vectorization  of  FASE  and  changes  to  the  input  and  output  routines.  Because  of  this, 
one  cannot  simply  use  the  ratio  of  total  run  times  to  calculate  differences.  Thus  the  ratios  given 
in  Table  2  must  be  combined  with  the  total  times  given  in  Table  3.  The  amount  of  time  saved  is 
given  by: 

0  SH  oL 
S=-r+-r 

T  T  Eq.  4 

where  d  and  s  are  given  as  "Savings"  in  Table  2,  H  and  L  are  the  total  time  spent  in  the  HERAC 
and  LBLF4  convolution  routines,  and  T  is  the  total  model  time  (FASE  2001)  from  Table  3. 

An  evaluation  of  the  amount  of  time  saved  by  eliminating  the  interpolation  of  the  coefficients 
(FASE  201  vs.  FASCODE),  by  using  fewer  points  to  define  the  shape  function  (FASE  2001  vs. 
FASE  201),  and  by  both  changes  (FASE  2001  vs.  FASCODE)  is  shown  in  Table  4.  The  total 
savings  should  be  the  sum  of  the  other  two  rows,  except  for  differences  between  FASCODE  and 
FASE  auxiliary  routines  (e.g.  BUFIN,  BUFOUT,  code  vectorization,  etc.).  The  time  saved  is 
relative  to  the  FASE  2001  times  given  in  Table  3.  The  "total”  values  for  S  in  Table  4  are 


somewhat  misleading  due  to  the  negative  values  for  the  time  "saved"  in  LBLF4.  This  is  more  of 
a  numeric  artifact  than  anything  else  since  the  individual  layer  times  for  LBLF4  are  either  very 
small  or  very  large,  but  are  always  very  similar.  (The  effect  of  the  interpolation  is  very  small 
relative  to  the  total  time;  the  number  of  points  to  describe  the  lineshape  is  irrelevant  in  this 
function).  While  Table  2  showed  relatively  large  savings  in  the  HIRAC  convolution  by  using 
fewer  points  or  by  eliminating  the  interpolation,  the  overall  effect  on  the  model  is  quite  small. 
This  is  true  particularly  in  Case  2  where  the  convolution  routines  are  a  much  smaller  percentage 
of  the  total  execution  time  than  for  the  other  cases  (~6%  for  Case  2  vs.  15-30%  for  Case  1  and 
Case  3). 


2.1. 1.2.  Voigt  Function  Numerical  Differences 

A  single-layer  case  (from  22  -  25  km)  was  used  to  examine  the  effects  of  additional  points  to 
define  the  lineshape  function  and  the  coefficient  interpolation.  The  results  of  this  test,  over  a 
narrow  portion  of  the  entire  spectral  region,  are  shown  in  Figure  1,  Figure  2,  and  Figure  3.  Figure 
1  shows  the  effect  of  interpolating  the  coefficients,  and  is  relatively  smooth  compared  to  the 
effect  of  increasing  the  number  of  points  used  to  describe  the  lineshape  function  (Figure  2).  The 
result  of  both  of  these  changes  is  shown  in  Figure  3,  where  the  overall  effect  produces 
differences  of  the  order  of  ±0.5%. 


Figure  1  FASE(201)  and  FASCODE  Differences 
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Figure  2  FASE(201)  and  FASE(2001)  Differences 
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Figure  3  FASE(2001)  and  FASCODE  Differences 
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It  is  clear  from  the  above  discussions  that  only  a  minimal  amount  of  time  can  be  saved  by 
reverting  to  the  FASCODE  lineshape  formulation.  Further,  while  the  optical  depth  differences 
between  the  "new"  FASE  algorithm  and  FASCODE  are  of  the  order  of  ±0.5%,  the  increased 
accuracy  is  required  for  many  remote-sensing  problems.  Thus  the  new  algorithm  does  not 
increase  the  execution  time  sufficiently  to  warrant  changes  in  this  algorithm. 

2.1.2.  Vectorization 

In  the  conversion  of  FASCODE  to  LBLRTM,  many  of  the  “do-loops”  were  restructured  to 
eliminate  “if-tests”  from  within  the  loops.  This  allows  the  FORTRAN  compiler  on  vector 
machines  to  make  optimal  use  of  the  machine  architecture,  thereby  decreasing  the  computation 
time.  A  side  benefit  is  that  the  code  itself  tends  to  be  easier  to  read,  modify,  and  debug. 

2.1.3.  Parameterization 

In  LBLRTM,  a  number  of  commonly  used  arrays  were  converted  from  “hard-wired” 
dimensions  to  dimensions  that  reside  in  parameter  statements.  Since  these  parameters  were 
given  the  same  names  in  different  modules,  it  is  relatively  easy  and  straightforward  to  change  the 
site  of  array  dimensions.  A  common  use  of  the  parameterization  is  to  increase  or  decrease  the 
number  of  atmospheric  layers  that  can  be  entered,  depending  on  the  size  of  the  input  profile  (e.g. 
radiosonde  data  can  contain  upwards  of  2000  points)  and  the  computer  memory. 

2.1.4.  Line  Rejection  Flag 

Line  rejection  refers  to  the  technique  of  rejecting  a  molecular  absorption  line  from  the 
calculation  if  its  strength  is  less  than  a  given  threshold.  Line  rejection  is  determined  separately 
for  each  layer  and  can  substantially  cut  the  computation  time  by  eliminating  a  large  number  of 
weak  lines  whose  effect  on  the  optical  depth  is  insignificant.  The  line  strength  threshold  is 
determined  by  the  requirement  that  the  optical  depth  at  the  line  center  be  less  than  some 
minimum  value  (the  default  is  0.0002.) 

A  line  rejection  flag  was  added  to  LBLRTM  to  signal  whether  or  not  a  particular  spectral 
line  was  rejected  for  use  in  a  layer.  This  binary  flag  is  stored  in  a  file,  which  can  be  read  in 
subsequent  runs  of  the  model  so  that  the  second  run  rejects  the  same  lines  as  the  first.  This 
feature  is  particularly  useful  when  computing  radiance  derivatives  with  respect  to  temperature 
using  a  finite-difference  scheme  since  a  line  could  be  rejected  in  the  reference  case  and  not  in  the 
test  case  after  perturbing  the  temperature  or  vice  versa.  In  such  cases  large  errors  would  be 
introduced  unless  the  line  was  consistently  rejected. 
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2.1.5.  Multiple  Scattering 

A  2-stream  approximation  for  multiple  scattering  was  developed  for  FASCODE  (Isaacs  et 
al.,  1987).  This  option  was  removed  from  LBLRTM  in  favor  of  an  option  to  create  the  files 
necessary  for  input  to  a  dedicated  multiple  scattering  code,  such  as  DISORT  or  CHARTS.  The 
rationale  was  that  many  applications  require  more  accuracy  than  that  provided  with  the  2-stream 
approximation,  and  that  a  stand-alone,  dedicated  multiple  scattering  code  is  able  to  provide  state- 
of-the-art  accuracy  with  minimized  computation  time. 

2.1.6.  LASER  Capability 

The  use  of  this  option  is  important  for  laser  calculations  because  of  the  FASCODE  4-function 
line  decomposition  algorithm  (see  discussion  of  spectral  grid  algorithm  in  Appendix  A).  In 
particular,  it  is  possible  that  a  specific  frequency  point  will  not  be  computed  at  all  altitudes,  but 
will  instead  be  determined  from  the  interpolation  of  adjacent  spectral  points.  The  LASER  option 
is  set  up  such  that  the  optical  depth  at  the  specific  LASER  frequency  is  computed  for  all 
atmospheric  layers,  thus  increasing  the  accuracy  of  the  calculation  for  this  specific  point.  This 
option  was  dropped  in  LBLRTM 

2.1 .7.  Non-LTE  Option 

Certain  molecular  species  can  exhibit  non-equilibrium  population  distributions  among 
vibrational  states  due  to  various  aspects  of  photochemistry  and  reduced  collision  rates.  This 
condition  is  known  as  non-local  thermodynamic  equilibrium,  or  NLTE,  and  mainly  occurs  at 
altitudes  above  40  km.  Molecules  that  are  ‘in  NLTE’  will  exhibit  thermal  radiative  emission  that 
can  not  be  described  by  the  Planck  function.  A  routine  for  NLTE  that  parallels  the  LTE  optical 
depth  routine  was  developed  for  FASCODE  (Ridgway  et  al.,  1982),  but  was  dropped  from 
LBLRTM. 

2.1.8.  Scanning  Function 

Many  applications  require  that  the  monochromatic  output  of  the  line-by-line  calculation  be 
convolved  (smoothed)  with  an  instrument  scanning  function.  FASCODE  performs  this  operation 
with  the  SCANFN  subroutine,  which  does  the  convolution  in  the  frequency  domain.  A  more 
efficient  and  accurate  method  is  to  use  Fourier  transform.  Let  S  be  the  monochromatic  spectrum, 
R  be  the  instrument  scanning  function  and  let  the  symbols  F  and  ★  represent  the  Fourier 
transform  and  convolution,  respectively.  A  fundamental  theorem  of  Fourier  transforms  states 
that: 


S  ★R  =  F,(F(S).F(R)) 


Eq.  5 
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Therefore,  the  convolution  can  be  performed  as  the  Fourier  transform  of  the  product  of  the 
transforms  of  the  spectrum  and  the  instrument  function. 

This  approach  has  been  implemented  in  the  subroutine  FFTSCAN  (Gallery  and  Clough, 
1992)  which  has  previously  been  available  only  as  a  stand-alone  program.  In  addition  to  being 
more  accurate  and  faster  than  SCANFN,  FFTSCAN  provides  for  a  variety  of  additional  scanning 
functions. 

2.1.9.  Plotting  Function 

The  plotting  function  was  originally  designed  for  creating  graphs  of  the  FASCODE 
output.  Given  the  proliferation  of  specialized  software  for  these  purposes,  the  main  utility  of  the 
plotting  function  is  to  create  a  2-column  ASCII  output  file. 

2.1.10.  Aerosol/Cloud  Attenuation 

The  cloud  and  aerosol  routines  found  in  FASCODE  and  LBLRTM  are  derived  from  the 
routines  found  in  LOWTRAN  (Anderson  et  ah,  1995).  Recent  improvements  to  these  routines 
for  MODTRAN  (Acharya  et  al.,  1993)  allow  for  multiple  cloud  layers  and  an  improved  interface 
for  user-supplied  information  (spectral  properties  and  vertical  profiles).  The  upgraded  routines 
are  discussed  below. 

2.1.11.  Spectral  Limit 

The  calculation  spectral  limit  was  changed  from  520  cm'1  in  FASCODE  to  1010  cm'1  in 
LBLRTM.  The  larger  limit  is  particularly  useful  in  the  mid-ER  where  instruments  typically  have 
a  large  spectral  range. 

2.1.12.  INPUT/OUTPUT  Routes 

The  BUFIN/BUFOUT  subroutines  control  the  unformatted  input  and  output.  The  coding  was 
changed  from  LBLRTM  to  increase  the  speed  of  the  I/O. 

2.1.13.  Geometry  Package 

All  three  codes  listed  in  Table  2-1  use  the  same  geometry  package  (Gallery  et  al.,  1983) 
which  includes  the  calculation  of  refractive  path.  This  package  was  found  to  have  some 
numerical  instabilities,  particularly  in  short  paths  that  are  nearly  tangent  to  the  atmosphere. 
These  errors  were  reduced  in  FASE  (see  Section  3.1). 
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2.2.  Numerical  Comparison 

A  comparison  of  FASCODE  and  LBLRTM  was  performed  by  AER  and  by  Dr.  Jinxue  Wang 
of  the  Phillips  Laboratory  to  ascertain  how  differences  in  the  model  physics  and  implementation 
(coding)  affect  the  final  model  output.  While  working  together  on  this  task,  these  groups  pursued 
different  strategies  for  the  comparison:  Dr.  Wang  chose  to  compare  the  FASCODE  and 
LBLRTM  radiances  to  ground-based  radiation  measurements  (ref),  while  the  AER  group  focused 
on  differences  in  the  calculated  optical  depths. 

There  are  several  differences  between  FASCODE  and  LBLRTM  that  should  be  examined 
before  discussing  the  simulation  studies.  First,  the  sub-Lorentzian  component  of  CO2  lines  is 
treated  differently.  This  gives  rise  to  numerical  differences  which  are  particularly  noticeable  near 
the  CO2  bandhead,  around  2380  cm*1,  and  one  would  expect  poor  agreement  between  the  models 
in  the  (approximate)  region  from  2375  -  2400  cm*1. 

The  implementation  of  line  coupling,  as  detailed  in  Section  4.3,  is  also  different  between 
FASCODE  and  LBLRTM.  The  coefficient  values  in  LBLRTM  were  computed  from  theory. 
These  coefficients  were  compared  with  measurements  and  scaled  by  30%  to  agree  with  the  data. 
This  scaling  is  not  inconsistent  with  the  calculations  because  it  is  within  the  uncertainty  of  the 
variables  in  the  calculations. 

Another  difference,  which  may  in  fact  be  the  most  important  difference  between  the  models, 
involves  the  monochromatic  spectral  resolution  (DV),  the  average  Lorentz  halfwidth  (ALFALO), 
and  the  maximum  halfwidth  (ALFMAX).  The  default  value  of  ALFALO  was  changed  from  0.08 
cm"1  in  FASCODE  to  0.04  cm'1  in  LBLRTM  because  0.04  cm'1  is  a  more  reasonable  number  for 
most  atmospheric  calculations.  The  value  of  ALFALO  impacts  the  resolution  (DV)  which  is 
defined  as  the  average  Voigt  halfwidth  (ALFV)  divided  by  the  number  of  samples  per  mean 
halfwidth: 


SAMPLE  Eq.  6 

where  ALFV  is  a  function  of  both  the  Doppler  and  Lorentz  halfwidths.  Changing  the  value  of 
ALFALO  by  a  factor  of  two  results  in  (almost)  a  factor  of  two  change  in  the  DV  and  increases 
the  number  of  spectral  points  in  the  LBLRTM  calculation  over  that  of  FASCODE.  Advances  in 
computer  memory  and  speed  since  the  original  development  of  FASCODE  means  that  increasing 
the  number  of  spectral  points  does  not  significantly  affect  the  calculation  time,  as  it  did  when 
0.08  cm'1  was  originally  chosen. 

Because  ALFALO  was  changed,  NBOUND  also  needed  to  be  changed  by  a  factor  of  two  in 
order  to  keep  ALFMAX  the  same: 
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ALFMAX  =  NBOUND* - * 


2  HWF3 

Eq.  7 

NBOUND= 

2  *(2  *  HWF3  *  SAMPLE+  0.0 1<_  FASCODE 

Eq.  8 

NBOUND= 

:  4  *(2  *  HWF3  *  SAMPLEf  0.0 1^^^^ 

Eq.  9 

While  ALFMAX  remains  the  same  after  these  changes,  there  is  still  a  factor  of  two  difference 
in  the  value  of  DV.  This  difference  between  the  two  models  can  produce  major  differences  in 
the  shape  of  the  spectrum:  DV  is  a  part  of  the  criteria  for  line  rejection  and  with  a  smaller  DV 
more  lines  will  remain  in  the  final  calculation.  However,  two  other  parameters,  DPTMIN  and 
DPTFAC,  also  influence  line  rejection.  Setting  these  parameters  equal  to  zero  will  ensure  that 
neither  program  rejects  any  lines. 

Another  aspect  of  the  difference  in  DV  is  that  the  monochromatic  spectra  from  FASCODE 
and  LBLRTM  will  be  written  to  TAPE  12  at  different  spectral  resolutions.  For  easy  comparison 
of  the  two  models  one  would  like  to  have  the  same  resolution.  One  hesitates  to  use  a  scanning 
function  to  change  the  DV  to  the  same  value  in  order  to  avoid  introducing  numerical  effects  into 
the  computed  spectrum.  Further,  it  is  important  that  the  user  not  use  the  DVSET  option  as  this 
will  fix  the  DV  during  the  calculations  and  can  introduce  distortions  to  the  shape  of  the  spectrum. 
Instead,  for  direct  comparison  of  LBLRTM  and  FASCODE  one  should  simply  interpolate  the 
final  LBLRTM  output  (which  is  at  higher  resolution)  to  the  resolution  of  FASCODE.  This  can 
be  done  using  the  INTRP  subroutine  or,  when  computing  optical  depths,  LBLRTM  can  be 
instructed  to  interpolate  the  output  to  a  specified  DV  without  changing  the  DV  used  in  the 
calculation  (via  the  subroutine  PNLINT  and  the  control  flag  DVOUT). 

Thus,  in  order  to  make  meaningful  comparisons  of  the  two  models  one  should  (a)  use  the 
default  values  for  ALFALO,  and  (b)  set  DPTMIN  and  DPTFAC  to  zero.  This  will  result  in  the 
DV  being  different  by  a  factor  of  two,  but  no  lines  will  be  rejected.  The  LBLRTM  spectra  can 
then  be  degraded  to  the  FASCODE  resolution  for  easy  comparison. 

A  minor  change  in  the  default  values  of  several  parameters  was  also  noticed.  The  default 
values  for  AVTRAT,  TDIFF1,  and  TDIFF2  were  changed  from  2.0,  8.0,  and  12.0  in  FASCODE, 
to  1.5,  5.0,  and  8.0  in  LBLRTM.  These  parameters  control  the  maximum  Voigt  width  ratio 
across  a  layer,  and  the  maximum  allowable  difference  in  temperature  across  the  layers. 
Decreasing  the  sizes  of  these  parameters  causes  an  increase  in  the  number  of  layers  used  in  the 
calculation  of  the  atmospheric  path.  This  can  provide  a  more  accurate  description  of  the 
atmosphere  at  the  expense  of  increased  computation  time. 

In  order  to  remove  extraneous  features  of  the  models,  which  can  complicate  the 
determination  of  differences  between  FASCODE  and  LBLRTM,  the  comparisons  were  based 


15 


solely  on  the  computation  of  transmission  for  a  single  atmospheric  layer.  This  was  accomplished 
for  three  different  layers  selected  to  represent  regions  of  predominately  Lorentz  broadening, 
Doppler  broadening,  and  an  intermediate  case.  Further,  only  the  first  seven  molecular  species 
from  the  HTTRAN  database  were  used,  and  the  column  density  for  each  of  the  molecules  was 
put  into  the  TAPE5  file  to  eliminate  differences  in  layering  and  the  calculation  of  molecular 
amounts.  Cross-sectional  and  continuum  data  were  not  included;  and  the  value  for  ALFALO 
was  set  to  the  default  value,  even  in  line-coupling  regions.  These  cases  correspond  to  conditions 
defined  by  the  US  Standard  Atmosphere.  Specific  characteristics  of  the  cases  are  summarized  in 
Table  5. 


Table  5  Parameters  for  FASCODE/LBLRTM  Comparisons 


Case  1 

Case  2 

Case  3 

Pressure  (mb) 

972.2107 

130.4209 

0.1251 

Temperature  (K) 

285.9446 

216.7000 

235.9188 

H20  (cm-2) 

1.231E+22 

6.600E+18 

3.810E+15 

C02  (cm-2) 

5.696E+20 

3.540E+20 

!  2.928E+17 

03  (cm-2) 

4.749E+16 

6.434E+17 

6.719E+14 

N20  (cm-2) 

5.524E+17 

3.176E+17 

1.424E+12 

CO  (cm-2) 

2.559E+17 

4.745E+16 

1.487E+14 

CH4  (cm-2) 

2.934E+18 

1.730E+18 

1.331E+14 

02  (cm-2) 

3.608E+23 

2.242E+23 

1.854E+20 

All  Others  (cm-2) 

1.349E+24 

8.514E+23 

7.281E+20 

Layer  Location  (km) 

0.0  -  0.7 

13.4-15.9 

62.9  -  65.3 

Approx.  Spectral  DV  (cm1) 

2.0E-02 

3.0E-03 

4.0E-04 

Because  of  various  array-sizes  in  FASCODE,  the  maximum  spectral  region  over  which  a 
calculation  can  be  done  is  520  cm'1.  Hence,  five  separate  runs  of  500  cm'1  regions  were  required 
in  order  to  cover  the  entire  spectral  range  from  500  -  3000  cm'1.  The  results  of  selected  cases  in 
this  spectral  region  are  shown  in  Part  (a)  of  Figure  4  through  Figure  6  as  the  difference  in 
transmission  between  FASCODE  and  LBLRTM.  For  information.  Part  (b)  of  these  figures 
illustrates  the  absolute  transmission  for  the  layer  computed  by  FASCODE  and  is  useful  in 
determining  optically  thick  or  thin  regions  of  the  spectrum  to  better  gauge  the  magnitude  of  the 
errors. 
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Transmission  Percent  Difference  in  Transmission 


Case  2 


Figure  5  Comparison  of  FASCODE  and  LBLRTM  Transmission  for  CASE2:  p=130.4209, 

T=216.7 
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In  general,  the  results  indicate  that  FASCODE  and  LBLRTM  agree  to  within  a  few  percent  in 
transmission,  which  corresponds  to  less  than  0.1%  in  optical  depth.  The  largest  differences 
among  Figure  4,  Figure  5,  and  Figure  6  occur  at  the  center  of  absorption  features.  In  order  to 
understand  the  cause  of  these  differences,  spectral  lines  were  examined  individually.  In  doing  so, 
it  was  found  that  differences  in  optical  depth  might  be  caused  by  the  way  the  models  compute  the 
lineshape.  To  reduce  computation  time,  both  models  compute  the  Voigt  lineshape  using  an 
accelerated  convolution  algorithm  (Clough,  et  al.,  1981)  which  decomposes  the  line  shape  into 
four  functions.  This  algorithm,  described  in  Section  2.1.1,  was  changed  for  LBLRTM  increase 
the  accuracy  of  the  calculation  as  compared  to  measured  spectra  (Clough  and  Brown,  1994).  The 
transmission  differences  seen  in  this  study,  which  are  usually  less  than  1%  in  optical  depth,  can 
be  attributed  to  the  change  in  the  Voigt  function.  As  illustrated  in  the  next  section,  these 
differences  have  little  effect  on  the  calculation  of  radiance  or  transmittance. 

2.3.  Conclusions 

From  the  code  comparisons  and  simulation  results  discussed  above,  we  concluded  that  the 
best  place  to  begin  the  FASE  model  would  be  from  LBLRTM,  since  it  already  has  the  most  of 
the  features  listed  in  Table  1.  The  main  features  it  lacks  are  the  non-LTE  capability  and  the 
internal  multiple  scattering  option.  However,  it  is  easier  to  put  these  back  into  the  model  than  to 
reintroduce  the  other  features  back  into  FASCODE.  Therefore  the  baseline  FASE  model  has  all 
the  LBLRTM  features  listed  in  Table  1 .  Further  modifications  were  made  to  FASE  to  upgrade 
existing  features  and  to  introduce  new  ones.  These  improvements  are  discussed  in  the  following 
sections. 
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3.  Upgrades  to  Existing  Features 


This  section  discusses  upgrades  to  features  already  existing  in  the  baseline  version  of 
FASE.  Some  of  these  upgrades  fix  errors  in  existing  code,  while  others  are  extensions  or 
improvements. 


3.1.  Geometry  Routines 

This  section  describes  changes  that  were  made  to  the  line-of-sight  geometry  routines.  As  an 
introduction  we  provide  a  brief  overview  of  the  input  parameters  and  how  the  package  computes 
the  final  parameters.  We  then  describe  problems  that  are  encountered  using  the  old  routines  and 
how  these  problems  have  been  solved  for  FASE.  This  work  closely  follows  changes,  which  were 
made  to  the  MODTRAN  (Anderson  et  al.  1994)  geometry  package,  as  detailed  in  the  report  of 
Acharya  et  al.  (1993).  Because  MODTRAN  uses  the  same  geometry  package  as  FASCODE,  and 
changes  designed  for  one  code  can  be  readily  ported  to  the  other  code. 

3.1.1.  Input/Output  Module 

The  TAPE5  input  file  will  accept  various  combinations  of  six  parameters  as  well  as  a  flag  for 
indicating  the  type  of  path.  These  parameters  are  listed  in  Table  6.  The  path  described  by  setting 
ITYPE  =  1  has  not  been  found  to  have  problems  and  will  not  be  discussed  here.  Paths  described 
by  ITYPE  =  3  are  a  special  case  of  the  ITYPE  =  2  path  and  thus  will  be  covered  by  changes  in 
the  geometry. 

Table  6  TAPE5  Input  Parameters 


Parameter 

Definition 

HI 

Observer  Altitude 

H2 

Source  Altitude 

ANGLE 

Zenith  Angle  at  H 1 

BETA 

Earth-centered  Angle  (HI  to  H2) 

RANGE 

Line-of-Sight  Distance  (H 1  to  H2) 

ITYPE  =  1 

Horizontal  Homogeneous  Path 

ITYPE  =  2 

Path  from  HI  to  H2 

ITYPE  =  3 

Path  from  H 1  to  Space 

For  the  ITYPE  =  2  specification  there  are  four  possible  combinations  of  input  parameters. 
These  are  listed  in  Table  7.  The  ultimate  goal  of  the  geometry  package  is  to  convert  the  input 
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parameters  into  the  starting  point  (HI,  which  is  always  known),  the  ending  point  (H2),  the 
viewing  angle,  and  the  minimum  point  along  the  path.  Thus,  all  paths  are  eventually  converted 
to  the  Case  2A  specification:  Case  2B  is  converted  to  2A  by  determining  H2;  Case  2C  is 
converted  to  2D  by  computing  the  earth-centered  angle  BETA;  and  Case  2D  is  converted  to  2A 
by  determining  the  angle. 

Table  7  ITYPE=2  Input  Parameters 


Case  Designation 

Input  Parameters 

2A 

HI,  H2,  ANGLE 

2B 

HI,  ANGLE,  RANGE  . 

2C 

HI,  H2,  RANGE 

2D 

HI,  H2,  BETA 

3.1.2.  Identification  of  Existing  Deficiencies 

Differences  between  the  path  specified  by  the  input  parameters  and  the  computed  line-of- 
sight  path  are  indicative  of  numerical  problems  within  the  geometry  routines.  As  indicated  in  the 
report  of  Acharya  et  al.  (1993)  the  inaccuracies  can  be  attributed  to  the  following: 

1 .  Numerical  precision  problems 

2.  Calculation  of  H2  without  including  refraction  (Case  2B) 

3.  Convergence  problems  for  calculation  of  BETA  (Case  2C  &  2D) 

4.  Short  slant  paths 

The  solutions  to  these  problems  are  discussed  separately  in  the  following  sections.  In 
addition  to  these  changes,  it  should  be  noted  that  an  attempt  was  made  to  make  code  easier  to 
read  by  eliminating  a  number  of  “GOTO”  statements. 

3.1.3.  Resolution  of  Deficiencies 
3.1 .3.1.  Numerical  Precision 

Many  of  the  equations  used  in  the  computation  of  the  geometry  involve  Re,  the  radius  of  the 
earth.  This  has  the  potential  for  introducing  numerical  errors  since  Re  is  usually  three  orders  of 
magnitude  larger  than  the  values  of  HI  and  H2.  In  order  to  increase  the  precision  of  these 
calculations  without  using  double  precision  variables,  equations  that  include  Re  were  re-written 
in  more  appropriate  forms.  For  example,  the  expression: 

R22+R2-R;  Eq.  10 
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can  be  reformulated  as 


H;-H,2  +  R2+2Re(H2-H,) 


Eq.  1 1 


since 


R,  =Re  +  H.  Eq.  12 

It  was  confirmed  that  the  numerical  result  of  Eq.  1 1  is  the  same  as  using  double  precision 
with  Eq  12.  Another  example  of  algebraic  reformulation  involves  trigonometric  identities: 


Rj+R^R^cosp 


Eq.  13 


becomes 


(H,  -H2)2  +4R1R,sin2 


(&) 

UJ 


Eq.  14 


3.1 .3.2.  Addition  of  Refraction  for  Calculation  of  Source  Altitude 

In  the  original  geometry  routine  H2  was  computed  from  HI,  ANGLE,  and  RANGE  by 
assuming  straight-line  geometry.  This  introduces  errors  that  are  especially  large  when  the  zenith 
angle  is  close  to  90°  and  the  range  is  large.  To  solve  this  problem  a  set  of  routines  was  added 
that  include  refraction  in  the  determination  of  H2.  This  is  accomplished  through  an  iterative 
procedure,  which  computes  the  line-of-sight  path  length  and  compares  it  with  the  desired  range. 
The  routine  converges  to  the  correct  value  for  H2  when  the  computed  range  is  equal  to  the  length 
of  the  input  range.  As  a  test  of  the  accuracy  the  range  was  re-computed  given  HI,  H2,  and 
ANGLE.  The  results  obtained  before  and  after  these  changes  were  made  to  the  code  were 
compared  with  the  results  obtained  using  MODTRAN2,  and  are  discussed  below.  We  did  not 
adopt  the  approach  of  MODTRAN2  and  change  all  of  the  geometry  routines  into  double 
precision.  However,  in  order  to  converge  to  the  results  of  MODTRAN2  it  was  necessary  to 
change  several  subroutines  to  double  precision.  Further,  many  variables  within  other  subroutines 
were  also  changed  to  double  precision.  These  additional  subroutines  and  changes  are  outlined  in 
Table  8. 
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Table  $  New  and  Modified  Routines 


New  Routines: 

NEWH2 

RTBIS 

FNDPTH 

Double-Precision  Routines: 

FINDSHDP 

SCALHTDP 

ANDEXDP 

RADREFDP 

Changed  Routines: 

RFPATH 

(some  changes  to  double  precision) 

ALAYER 

FNDHMN 

3.1 .3.3.  Short  Slant  Paths 

The  report  of  Acharya  et  al.  (1993)  recommended  a  new  algorithm  (based  on  the  Newton- 
Raphson  iteration  method)  for  the  calculation  of  BETA.  While  the  implementation  of  this 
MODTRAN2  routine  should  have  been  straightforward,  this  was  not  the  case.  Due  to  slight 
numerical  differences  between  MODTRAN  and  FASE  the  routine  would  not  converge  to  the 
correct  solution.  In  the  end  it  was  decided  that  the  old  "FDBETA"  routine  coupled  with  the 
algebraic  and  double-precision  modifications  discussed  in  Section  3. 1.3.1  was  sufficiently  better 
than  the  original  geometry  package.  Similarly,  no  modifications  were  made  to  improve  the 
outputs  for  short  slant  paths. 

3.1 .3.4.  Refractive  Index  Profile 

In  the  course  of  validating  the  changes  to  the  geometry  it  was  noticed  that  MODTRAN2  and 
FASE  use  different  formulations  for  the  index  of  refraction.  Further  investigation  showed  that 
both  formulations  are  based  on  the  paper  by  Edlen  (1966).  The  formulation  in  FASE  (and 
FASCODE  and  LBLRTM)  was  originally  coded  for  LOWTRAN  and  is  presented  in  the  report 
for  LOWTRAN5  (Kneizys  et  al.,  1980).  This  formulation  is  the  simplification  of  Edlen’s 
expression.  Apparently  the  simplification  was  dropped  in  LOWTRAN6  in  favor  of  the  actual 
expression  (Kneizys  et  al.,  1983),  and  thus  was  employed  in  MODTRAN.  While  the  profiles 
exhibit  only  slight  differences  these  are  sufficient  to  cause  small  deviations  in  the  final  line-of- 
sight  path,  especially  for  cases  where  refraction  is  significant.  Since  the  LOWTRAN5 
expression  is  a  simplification  to  the  exact  expression  we  have  dropped  it  from  the  coding  in  favor 
of  the  exact  expression  (this  change  was  also  made  in  the  subroutine  XAMNTS).  We  do  not 
anticipate  that  this  change  will  require  the  revalidation  of  FASE  with  other  models  or 
measurements. 
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3.1.4.  Validation  of  Geometry  module 

In  order  to  gauge  the  impact  of  the  coding  changes,  comparisons  were  made  between  the 
"new"  FASE  geometry,  the  "old"  FASE  geometry,  and  the  MODTRAN2  geometry.  The  results 
of  these  comparisons  are  shown  in  Table  9  for  Case  2A,  Table  10  for  Case  2B,  and  Table  1 1  for 
Case  2C/2D. 

Case  2A  modifications  consisted  of  the  reformulation  of  equations  and  the  changing  of  some 
variables  to  double  precision.  Therefore,  one  would  not  expect  large  differences  between  the  old 
and  new  geometry.  Table  9  illustrates  changes  to  the  value  of  phi,  the  zenith  angle  at  H2  (used  in 
the  calculation  of  the  refracted  path),  as  a  result  of  the  new  coding.  The  coding  changes  had  little 
impact  on  the  final  results. 


Table  9  Case  2A  Model  Comparisons  for  Geometry 


HI 

H2 

ANGLE 

PHI 

(old) 

PHI 

(new) 

HMIN 

(old) 

HMIN 

(new) 

20 

100 

90.1 

99.0120 

99.0120 

19.9902 

19.9901 

90.5 

99.0252 

99.0252 

19.7518 

19.7521 

91 

99.0663 

99.0663 

19.0059 

19.0057 

60 

100 

90 

96.3737 

96.3737 

60.0000 

60.0000 

90.1 

96.3745 

96.3745 

59.9902 

59.9900 

90.5 

96.3932 

96.3932 

59.7549 

59.7551 

91 

96.4514 

96.4514 

59.0205 

59.0206 

95 

98.0945 

96.0945 

35.5166 

35.5167 

96 

98.7450 

98.7450 

24.7101 

24.7101 

Larger  differences  between  the  old  and  new  geometry  are  apparent  for  Case  2B,  shown  in 
Table  10.  The  largest  errors  in  the  old  geometry  occur  when  the  refraction  is  largest,  at  viewing 
angles  close  to  90°.  Also,  the  new  geometry  will  give  results  for  the  case  of  a  90°  viewing  angle, 
while  MODTRAN  returns  an  error  message. 
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Table  10  Case  2B  Model  Comparisons  for  Geometry  (N,0  =  new, old  FASE,  M  = 
MODTRAN;  for  PHI  =  90°  MODTRAN2  is  unable  to  find  the  tangent  point) 


Model  HI 

ANGLE 

RANGE 

H2 

PHI 

RANGE 

92 

10 

4.65794 

88.0803 

10.01368 

4.65794 

88.0803 

10.00394 

4.65869 

88.0794 

9.97758 

N  5 

92 

50 

3.42779 

88.3960 

50.0190 

M 

3.42779 

88.3959 

50.0214 

0 

3.45068 

88.3887 

49.1827 

N  5 

92 

100 

2.19749 

88.78511 

100.0631 

M 

2.19743 

88.7849 

100.04798 

0 

2.29346 

88.7495 

95.5356 

N  5 

92 

500 

4.46282 

91.8740 

503.416 

M 

4.49087 

91.8811 

504.337 

0 

7.15381 

92.4470 

574.756 

N  5 

89 

10 

5.18148 

91.0812 

10.00415 

M 

5.18149 

91.0812 

10.00076 

O 

5.18213 

91.0799 

10.02902 

N  5 

89 

50 

6.04695 

91.4015 

49.9906 

M 

6.04695 

91.4015 

49.9993 

O 

6.06836 

91.4078 

50.8677 

N  5 

89 

100 

7.44533 

91.8068 

99.9934 

M 

7.44527 

91.8068 

99.9978 

O 

7.52881 

91.8271 

102.6277 

N  5 

89 

500 

31.8327 

95.2424 

499.9550 

M 

31.8336 

95.2422 

499.9918 

O 

33.2676 

95.3797 

515.4590 

N  5 

90 

10 

5.00715 

90.0947 

10.16462 

0 

5.00781 

90.0862 

10.60430 

N  5 

90 

50 

5.17290 

90.3990 

49.9615 

0 

5.19580 

90.4229 

53.0766 

N  5 

90 

100 

5.62983 

90.7967 

99.9713 

0 

5.78418 

90.8484 

106.150 

N  5 

90 

500 

22.7348 

94.1480 

499.950 

0 

24.5737 

94.3697 

523.980 
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The  results  for  Case  2C/2D,  shown  in  T able  1 1 ,  also  represent  an  improvement  over  the  old 
routine,  which  would  not  converge  at  all  for  the  case  shown. 

Table  11  Case  2C/2D  Model  Comparisons  for  Geometry  (Hl=H2=5km) 


Input  Range 

Output  Range 

FASE 

Output  Range 

MODTRAN 

Output  Range 

MODTRAN2 

300 

300.044 

300.13 

300.01 

200 

199.869 

199.96 

200.01 

100 

100.3059 

100.52 

100.01 

50 

49.6897 

50.51 

50.02 

20 

19.8501 

- 

20.01 

10 

9.35703 

9.20 

10.02 

8 

8.10367 

7.51 

8.01 

6 

- 

7.51 

6.01 

4.7 

- 

- 

4.72 

2.01 

0.0 

5.31 

2.00 

Table  12  summarizes  the  geometry  test  results  for  FASE,  MODTRAN2,  and  FASCOD3P. 
For  this  test  HI  =  5.0  km.  Range  =  10  km,  and  H2  was  varied  as  indicated  in  the  table.  Given 
HI,  H2,  and  Range  the  geometry  package  will  first  calculate  the  earth-centered  angle  ("PHI")  and 
then  the  H1-H2  viewing  angle  ("ANGLE”).  The  numbers  tabulated  are  the  range  computed  after 
finding  the  viewing  angle  (and  should  be  equal  to  the  input  range  of  10  km).  Except  for  Test  2, 
the  new  FASE  geometry  is  much  better  than  the  old  FASCODE  geometry.  However,  as  noted  in 
the  Version  1.0  report,  there  are  still  some  differences  relative  to  MODTRAN. 

Further  examination  of  these  test  cases  indicates  that  the  differences  between  FASE  and 
MODTRAN  are  the  result  of  slight  numerical  differences. 

Table  13  shows  the  values  of  PHI  (the  zenith  angle  at  H2  used  in  the  calculation  of  the 
refracted  path),  ANGLE  (the  viewing  angle  at  HI),  and  HMIN  (the  computed  tangent  point  along 
the  path  from  HI  to  H2).  "RANGE  DIFF”  refers  to  the  difference  in  the  range  computed  with 
FASE  and  MODTRAN2  as  given  in  Table  12;  "PHI"  is  the  earth-centered  angle;  "ANGLE"  is 
the  H1-H2  viewing  angle;  and  "HMIN"  is  the  computed  tangent  point  of  the  refracted  path.  The 
results  of  Test  3  clearly  show  that  just  a  slight  change  in  the  calculation  of  the  viewing  angle  will 
affect  the  final  calculation  of  the  range. 
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Table  12  Case  2C  Geometry  Test  Results 


H2 

FASE 

FASCOD3P 

MODTRAN2 

6.0 

10.001 

10.008 

5.5 

10.007 

10.005 

5.25 

10.012 

10.024 

10.014 

5.0 

9.357 

11.837 

10.014 

5.01 

10.335 

no  convergence 

10.009 

5.007 

10.234 

9.547 

-  I 

7 

5.001 

no  convergence 

no  convergence 

'  | 

5.0015 

10.219 

no  convergence 

10.013 

Table  13  Values  Computed  to  Define  Path  Geometry 


TEST 

RANGE  DIFF 

FASE 

MODTRAN2 

3 

0.002 

PHI 

91.470 

91.471 

ANGLE 

88.609 

88.609 

HMIN 

5.000 

5.000 

5 

0.326 

PHI 

90.096 

90.097 

ANGLE 

89.986 

89.983 

HMIN 

5.000 

5.000 

8 

0.206 

PHI 

ANGLE 

90.037 

90.031 

HMIN 

4.999 

4.999 

In  summary,  the  results  of  the  validation  study  indicate  that  the  new  geometry  package  is  a 
significant  improvement  over  the  existing  routines  and  solves  many  of  the  problems  that  have 
been  encountered.  Further  changes  to  the  geometry  to  improve  Case  2C/2D,  especially  for  short 
input  ranges,  will  be  considered  for  future  versions  of  FASE  as  time  permits. 

3.2.  Chappuis  And  Wulf  Ozone  Bands 

The  coefficients  for  the  ozone  Chappuis  band  were  updated  to  the  values  found  in 
MODTRAN3  (Shettle  and  Anderson,  1995).  These  coefficients  include  both  the  Chappuis  and 
Wulf  absorption  and  cover  the  range  from  9170  to  24565  cm*l.  The  coefficients  at  each 
wavenumber  are  now  given  for  a  quadratic  polynomial  that  is  a  function  temperature: 
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C(u)  =  X(d)  + Y('U)A  +  Z(d)A2 


Eq.  15 


where 


A  =  T  -273.15  Eq.  16 

and  T  is  the  average  temperature.  The  coefficients  are  specified  at  intervals  of  5  cm-'. 

The  original  coefficients  were  valid  over  the  temperature  range  220  -  298  K.  This  was 
extended  to  180  -  310  K  to  cover  the  range  of  temperatures  typically  encountered  in  atmospheric 
applications.  Further,  an  extrapolation  was  done  at  the  edges  of  the  data  to  eliminate  edge  effects 
caused  by  a  sudden  decrease  of  the  coefficients  to  zero  value. 

3.3.  Non-Local  Thermodynamic  Equilibrium  (NLTE) 

One  of  the  larger  projects  in  converting  LBLRTM  into  FASE  involved  the  conversion  of  the 
FASCODE  NLTE  routines  into  a  format  compatible  with  the  LTE  portion  of  the  code  (panel 
structure,  common  blocks,  array  dimensions,  etc.).  Prior  to  making  this  conversion,  several 
differences  between  the  LTE  and  NLTE  routines  of  FASCODE  were  identified  and  eliminated 
(by  Jim  Chetwynd  of  the  Air  Force  Phillips  Laboratory)  so  that  a  calculation  using  only  LTE 
spectral  lines  will  achieve  the  same  answer  regardless  of  whether  the  LTE  or  NLTE  routines 
were  used.  Once  these  changes  were  made  to  FASCODE,  the  sections  of  code  relevant  to  NLTE 
calculations  were  identified.  These  sections  were  then  added  to  the  FASE  subroutines  to  create  a 
parallel  version  of  HIRAC  for  NLTE  calculations. 

3.4.  Black-Body  Approximation 

At  the  October  1994  ARM  Science  Team  meeting  it  was  shown  that  there  were  problems  in 
the  LBLRTM  output  when  doing  calculations  with  aerosols,  as  evidenced  by  a  "step  function" 
radiance  profile  rather  than  a  smooth  curve  (Anderson,  1994).  We  worked  directly  with  the 
LBLRTM  team  at  AER  to  identify  the  cause  of  these  problems  and  to  re-code  the  necessary 
sections  of  the  model.  It  was  found  that  the  step-function  was  the  result  of  improper 
interpolations  of  the  blackbody  function  as  a  part  of  the  Pade  approximation.  This  error  was 
related  to  another  problem  in  which  LBLRTM  did  not  properly  include  the  surface  emissivity. 
We  have  corrected  these  sections  of  the  code  in  both  FASE  and  LBLRTM. 

3.5.  Emissivity/Reflectivity 

Previous  versions  of  FASE,  FASCODE,  and  LBLRTM  give  the  option  of  specifying  the 
boundary  emissivity  and/or  reflectivity.  This  is  accomplished  through  the  use  of  quadratic 
coefficients  to  specify  the  spectral  signature  of  the  boundary.  However,  a  quadratic  formulation 
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is  not  sufficient  for  all  simulation  scenarios  since  the  emissivity  or  reflectivity  can  have  a  non- 
uniform  spectral  distribution  (e.g.,  the  sea-surface  emissivity).  To  rectify  this  problem  we  have 
included  the  option  of  putting  the  emissivity  and  reflectivity  coefficients  in  data  files.  This  is 
similar  to  an  option  that  was  added  to  MODTRAN.  To  use  this  option,  set  the  first  coefficient 
(SREMIS(l)  and/or  SRREFL(l))  equal  to  -99  and  the  other  two  coefficients  equal  to  zero.  This 
will  signal  FASE  to  look  for  the  files  "EMISIV.IN"  and/or  "REFLEC.IN". 

3.5.1 .  Input  File  Format 

The  format  for  the  emissivity  and  reflectivity  files  is  the  same.  The  first  line  contains  header 
information  that  will  be  written  to  the  TAPE6  output  file.  The  second  line  gives  the  starting  and 
ending  wavenumbers  for  the  file  as  well  as  the  number  of  points  in  the  file.  The  user  may  set  the 
starting  wavenumber  to  be  less  than  zero  to  signify  that  wavenumber-coefficient  pairs  are 
specified,  otherwise  it  is  assumed  that  the  coefficients  are  spaced  equally  between  the  starting 
and  ending  wavenumbers.  The  remaining  lines  in  the  data  file  contain  the  coefficients  or 
wavenumber-coefficient  pairs.  These  are  listed  with  six  values  per  line  (or  three  pairs  if  the 
wavenumber  option  has  been  selected).  This  file  structure  is  summarized  below: 

line  1:  80  character  header  (for  user  information  only  -  written  to  TAPE6) 

line  2:  VI,  V2,  NPTS  (E13.6,  E13.6, 14) 

VI  =  starting  wavenumber  of  data 

V2  =  ending  wavenumber  of  data 

NPTS  =  number  of  data  points  in  file  (max  =  2000) 

line  3  -  n:  coefficients  ( 6(E11.4,2X) ) 

coefficients  are  listed  with  six  values  per  line,  with  equal  wavenumber  spacing 
between  the  coefficients 

if  V 1  <  0  in  line  2,  the  wavenumber-coefficient  pairs  are  expected  and  there  are 
three  pairs  per  line 

3.6.  Hartley-Huggins  Continuum  /  Herzberg  Continuum 

An  error  in  the  array  dimensions  for  the  Hartley-Huggins  (O3)  continuum  absorption  data 
statements  was  identified  and  corrected. 

The  pressure-dependence  correction  routine  for  the  Herzberg  (O2)  continuum  cross-sections 
has  been  modified  to  reflect  a  more  recent  formulation  (Anderson  et  al.,  1990).  The  pressure 
dependence  proportionality  constant  was  changed  from  1.72E-03  (at  760  torr  and  273.15  K)  to 
1.81E-03  (at  760  torr  and  293.15  K).  The  net  result,  after  converting  to  273.15  K  and  including 
the  O2  and  N2  mixing  ratios,  is  a  change  from  0.73  to  0.83  in  the  factor  modifying  the 

pressure/temperature  scaling  term. 
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3.7.  Absorption  Cross-Sections 

A  number  of  "go  to"  statements  in  the  subroutine  "XSECTM"  were  changed  to  avoid  an  error 
whereby  all  remaining  cross-section  files  would  be  skipped  if  it  was  found  that  the  current  file 
does  not  have  any  lines  in  the  spectral  region  of  interest.  It  was  also  discovered  that  some  array 
indices  were  incorrectly  coded  after  allowing  for  the  option  to  input  cross-sections  with  pressure 
units  of  Torr. 

3.8.  Layer  Input/Output  Units 

One  of  the  options  available  is  for  the  user  to  specify  the  layering  scheme  and  provide  the 
column  amounts  within  each  layer  for  the  molecules  of  interest.  We  have  extended  the  allowed 
input  units  to  include  "atm-cm".  If  this  option  is  selected  the  input  amounts  are  converted 
directly  to  molecules  cm' 2;  without  this  option  the  input  is  assumed  to  be  molecules  cm'2  unless 
the  value  is  less  than  1.0  in  which  case  it  is  assumed  to  be  mixing  ratio.  With  the  mixing  ratio 
option  the  column  amount  for  the  broadening  gas  must  still  be  input  in  molecules  cm~2,  while  for 
the  atm-cm  input  all  of  the  amounts  must  be  in  atm-cm. 

Using  the  IPUNCH  =  1  option  (TAPE5,  record  3.1),  layer  information  written  to  TAPE6  will 
now  include  the  units  of  atm-cm. 

3.9.  Laser  Calculations 

As  noted  in  Section  2.1.6,  the  4-function  line  decomposition  employed  by  FASCODE  can 
introduce  errors  in  the  calculation  of  a  specific  frequency  point,  and  the  LASER  option  was 
devised  to  circumvent  this  for  calculations  over  a  very  narrow  spectral  band.  Since  the  laser 
option  has  been  recoded  for  use  in  FASE,  it  was  dropped  from  LBLRTM.  This  option  allows  the 
user  to  calculate  only  the  narrow  region  around  the  laser  frequency  and  print  the  results  at  the  end 
of  the  TAPE6  output  file.  This  option  is  signaled  if  V1=V2  in  the  TAPES  input  file. 
Calculations  with  the  laser  option  were  compared  to  calculations  over  a  wide  spectral  region  of 
radiance,  transmittance,  and  optical  thickness  and  were  found  to  be  in  exact  agreement. 
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4.  New  Features 


This  section  describes  new  features  added  to  FASE  that  were  not  contained  in  any  of  the 
previous  FAS  CODE  models  or  in  LBLRTM. 

4.1.  Schumann-Runge  Band  and  Continuum 

Information  about  the  02  Schumann-Runge  band  and  continuum  absorption  was  added  to  the 
continuum  module  of  FASE.  The  data  are  polynomial  coefficients  from  a  fit  devised  by 
Minschwaner  et  al.  (1992,  1995)  and  are  valid  from  49000  cm'1  to  57000  cm'1.  The  fit  was 
accomplished  on  a  0.5  cm'1  spectral  grid  using  cross-sections  obtained  from  a  line-by-line 
radiative  transfer  model  which  included  contributions  from  the  temperature  dependent 
Schumann-Runge  continuum. 

For  a  single  layer  the  absorption  cross-section  at  a  particular  wavenumber  may  be  written  in 
the  form 


S  =  Ax2  +Bx  +  C 


Eq.  17 


where  S  has  units  of  cm2  and  x  is  the  "modified  temperature",  which  is  related  to  the  average 
temperature  of  the  layer,  T: 


T-IOOY 

10  J  ' 


Eq.  18 


In  order  to  achieve  the  proper  dependence  on  temperature  there  are  three  sets  of  coefficients 
corresponding  to  three  temperature  regimes.  These  regimes  are  defined  as: 

130  K  •  T  •  190  K 
190  K  <  T  •  280  K 
280  K  <  T  •  500  K 


For  convenience  the  data  are  supplied  in  three  ASCII  data  files  (130-190.cf4,  190-280.cf4, 
and  280-500.cf4). 

The  cross-section  coefficients  are  tabulated  at  intervals  of  0.5  cm  .  For  three  coefficients  in 
three  temperature  regimes  this  corresponds  to  a  total  of  72000  values.  Rather  than  place  all  of 
these  values  in  "block  data"  statements,  a  program  was  written  to  convert  the  ASCII  data  files 
into  an  unformatted  data  file  (SCHRUN.DAT)  which  can  be  read  using  the  FASE  I/O  subroutine 
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"BUFIN".  Because  unformatted  data  is  treated  differently  with  different  computers,  FASE  is 
distributed  with  ASCII  data  files  and  a  conversion  program. 

The  unformatted  data  file  is  written  in  a  block  structure  consisting  of  a  header  and  the 
corresponding  data.  The  header  contains  the  starting  wavenumber  of  the  block,  the  spectral 
spacing  of  the  data,  the  number  of  points,  and  a  flag  signifying  the  temperature  regime.  Each 

block  of  data  covers  400  cm  and,  with  three  coefficients  for  each  spectral  point  and  a  resolution 

-1 

of  0.5  cm  ,  each  block  contains  2400  points. 

Within  the  continuum  module  the  program  chooses  the  correct  temperature  regime,  reads  the 
corresponding  coefficients,  and  computes  the  cross-section  as  in  Section  4.5.  The  result  is 
interpolated  to  the  appropriate  spectral  resolution  and  added  to  the  array  containing  the 
continuum  contribution  to  the  total  optical  properties. 

4.2.  Solar  Spectrum 

A  solar  spectrum  module  was  added  to  give  FASE  the  capability  of  computing  solar 
transmission  through  the  atmosphere  (which  is  particularly  useful  for  simulations  of  experiments 
employing  the  technique  of  solar  absorption  spectroscopy).  Ultimately  the  solar  module  could 
be  combined  with  a  multiple-scattering  routine  for  a  better  treatment  of  the  atmospheric  radiation 
field.  FASE  is  distributed  with  two  solar  spectrum  files  to  cover  the  spectral  range  from  50  - 
57490  cm'1. 

The  spectrum  was  supplied  in  SOL. RAD  by  Kurucz  (1994).  It  is  composed  of  three  solar 
irradiance  spectra  collected  by  Kurucz,  the  Air  Force  Phillips  Laboratory  Geophysics  Directorate 
(GL),  and  the  SUSIM  instrument  which  flew  on  the  Space  Lab  2  flight  of  the  Space  Shuttle  (Hall 
and  Anderson,  1991).  The  following  details  the  method  by  which  the  data  were  converted  to  the 
spectrum  supplied  with  FASE. 

The  Kurucz  spectrum  was  given  as  the  value  of  the  solar  irradiance  at  intervals  of  1  cm'1  from 
50  to  50000  cm  1  in  units  of  ergs  cm'2  s'1  (cm1)1.  The  irradiance  was  converted  to  the  standard 
unit  of  watts  cm"  (cm1)1,  and  the  spectrum  was  then  smoothed  using  a  running  three-point 
triangle,  to  make  its  excursions  more  similar  to  measured  spectra  in  the  ultraviolet.  It  should  also 
be  noted  that  toward  higher  wavenumbers  the  smoothing  begins  to  have  little  effect,  since  the 
spectrum  becomes  naturally  smooth. 

The  GL  spectrum  is  a  combination  of  1978  data  and  1983  data,  normalized  to  1983.  Given  in 

O 

0.02  A  steps  and  was  also  smoothed  by  a  3-point  triangle,  it  was  then  normalized  to  the  SUSIM 
data  from  the  Space  Lab  2  flight  using  a  41 -point  smoothing  of  the  ratio  of  the  two  spectra.  The 
data  were  then  converted  from  irradiance  units  of  photons  cm"  s'2  A'1  to  watts  cm'2  (cm1)  '*  and 
interpolated  to  give  values  at  equally  spaced  intervals  of  1  cm1. 
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The  Kurucz  spectrum  and  the  GL  spectrum  were  then  compared  in  the  range  of  32250-32350 
cm'1  and  a  splice-point  was  chosen  on  an  ascending  flank  in  the  range  from  32288  to  32292  cm'1, 
where  the  agreement  was  good  for  five  points  in  a  row.  The  composite  spectrum  is  the  smoothed 
Kurucz  spectrum  from  51  to  32290  cm1  and  the  smoothed  and  interpolated  GL  spectrum  from 
32291  to  49983  cm'1. 

The  second  spectrum,  SOL2.RAD,  covering  the  range  from  48560  to  57490  cm'1  is  from 
SUSIM  data.  It  has  the  same  units  and  file  structure  as  SOL.RAD. 

4.2.1.  File  Structure 

The  solar  spectrum  data  were  initially  supplied  in  a  2-column  ASCII  data  file.  This  file  is 
rather  large  (about  898  K)  and  the  first  step  was  to  reduce  it  to  a  more  manageable  size.  This  was 
accomplished  by  converting  it  to  a  block  format  (still  in  ASCII).  The  blocks  consist  of  a  header 
denoting  the  starting  wavenumber  of  the  block,  the  spectral  interval  (1  cm1),  and  the  number  of 
points  in  the  block  (2400  for  all  except  the  last,  which  contains  1933  points).  The  irradiance  data 
follows  each  header  and  is  arranged  with  six  values  per  line.  (The  very  last  line  has  one  real 
point  and  five  denoted  "-99.9"  to  fill  the  line).  Converting  the  data  to  this  format  reduces  the  file 
size  to  about  608  K. 

As  with  the  data  for  the  Schumann-Runge  band  (Section  4.1)  it  is  advantageous  to  read  data 
from  an  unformatted  data  file.  The  user  is  supplied  with  the  Kurucz  solar  spectrum  in  the  block- 
ASCn  format  (SOL.ASC).  A  conversion  program  is  supplied  that  will  convert  the  block-ASCH 
file  into  an  unformatted  file  (SOL.UNF)  with  the  same  structure  (this  program  ignores  the  last 
five  data  values,  "-99.9",  to  ensure  that  they  are  not  inadvertently  incorporated  into  FASE 
calculations).  The  unformatted  file  occupies  about  200  kilobytes  of  disk  space.  The  SOL2.UNF 
file  is  created  with  another  conversion  program  converting  block-data  statements  into  an 
unformatted  file. 

Supplying  the  solar  spectra  in  this  fashion  is  advantageous  to  the  user  community  since  each 
user  may  easily  use  any  solar  data  without  converting  it  to  a  block-data  statement  and 
recompiling  the  program,  thereby  increasing  the  flexibility  of  FASE. 

4.2.2.  Implementation 

Implementation  of  the  solar  spectrum  is  accomplished  with  Record  1.2  and  Record  1.2b  of 
the  TAPES  input  file.  The  user  has  the  choice  of  doing  a  ’normal’  calculation  and  then 
incorporating  the  solar  spectrum,  or  simply  using  previously  calculated  files.  In  order  to  do  the 
complete  calculation,  EMIT  (Record  1.2)  should  be  set  equal  to  2;  EMIT  should  be  3  to  skip 
directly  to  the  solar  routines.  With  either  case  the  filter,  scan,  and  plot  functions  occur  after  the 
solar  routines. 
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FASE  requires  that  Record  1.2b  follow  Record  1.2  if  EMIT  equals  2  or  3  (note  that  Record 
1.2  is  the  first  line  after  the  header  information).  Record  1.2b  contains  three  parameters 
(INSOLR,  IOTSOL,  and  ISOLOC)  which  provide  information  about  input  and  output  files,  and 
thus  the  type  of  merge  that  is  to  be  done.  If  INSOLR  equals  0  the  input  file  is  the 
radiance/transmittance  file  (TAPE  12),  while  if  it  equals  1  the  input  file  is  the  total  optical  depth 
file  (TAPE10).  With  IOTSOL  equal  to  0  the  merge  produces  only  the  attenuated  solar  radiation, 
while  IOTSOL  equal  to  1  produces  the  total  radiation  field  (the  transmitted  solar  radiation  plus 
the  thermal  emission).  The  output  is  always  written  to  TAPE1 1.  One  should  also  recognize  that 
if  the  total  radiation  field  is  desired,  INSOLR  must  be  equal  to  zero  to  provide  the  thermal 
emission  term.  The  ISOLOC  parameter  is  used  to  specify  the  location  of  the  solar  irradiance.  If 
ISOLOC  equals  0  (default),  the  solar  irradiance  is  taken  from  the  file  SOL.UNF  (or  SOL2.UNF), 
while  for  ISOLOC  equal  to  1  the  solar  irradiance  is  read  from  TAPE24.  The  TAPE24  file  is 
created  by  setting  ISOLOC=l,  as  discussed  in  Section  4.2.3. 

The  merging  of  the  solar  radiance  with  the  input  information  is  accomplished  as  follows:  If 
optical  depth  is  input,  the  transmission  is  computed;  the  solar  file  is  interpolated  to  the  spectral 
resolution  of  the  input  file  and  the  transmitted  solar  radiation  is  computed.  Finally  the  thermal 
component  is  added  (if  desired)  and  the  result  is  written  to  TAPE1 1. 

4.2.3.  Calculations  with  the  Solar  Irradiance  Option 

Calculations  utilizing  the  solar  irradiance  option  fall  into  two  categories.  The  first  is  direct 
transmission  where  the  observer  is  looking  at  the  sun,  either  with  a  short,  direct  slant  path  or 
through  a  limb  path.  The  other  type  of  calculation  is  that  of  sunlight  reflected  off  a  surface  into 
the  line-of-sight. 

Depending  on  the  geometry  of  the  problem,  calculations  with  a  surface  reflectance  can 
require  multiple  runs  of  FASE.  For  a  downlooking  calculation,  FASE  sets  the  top  of  the 
atmosphere  to  the  observer  altitude  (HI)  and  computes  both  the  upwelling  radiance  from  the 
surface  (H2)  to  the  observer  and  the  down  welling  radiance  from  HI  to  H2,  reflected  off  the 
surface  and  attenuated  back  to  HI.  Whether  or  not  this  calculation  is  sufficient  depends  on  the 
altitude  of  HI  and  H2  (both  relative  and  absolute)  and  the  spectral  band.  In  a  weakly  absorbing 
spectral  region,  the  downlooking  instrument  at  HI  may  'see'  radiation  that  was  emitted  from 
above  HI,  attenuated  down  to  the  reflective  surface  H2,  reflected  and  attenuated  back  to  HI. 
These  types  of  calculations  can  be  done  with  FASE,  but  require  varying  degrees  of  setup  by  the 
user. 

The  downlooking  case  from  HI  to  a  reflective  surface  H2  with  no  contribution  to  the 
radiance  from  the  atmosphere  above  HI  can  be  run  in  the  same  manner  as  the  direct  transmission 
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case:  set  IEMIT=2,  INSOLR=0,  IOTSOL=0  or  1,  and  ISOLOC=0.  One  should  note  that  the 
surface  is  treated  as  a  Lambertian  reflector  for  the  purposes  of  the  reflected  solar  irradiance. 

A  calculation  that  needs  to  include  radiation  emitted  from  above  HI  must  be  done  as  two 
runs  of  FASE.  For  purposes  of  this  discussion  the  observer  altitude  will  be  defined  as  Hobs,  the 
surface  altitude  as  Hsurf,  and  the  top  of  the  atmosphere  (typically  100  km)  as  TO  A.  The  first 
run,  "run  1”,  should  be  set  as  for  the  uplooking  direct  transmission  case:  Hl=Hobs,  H2=TOA, 
IEMIT=2,  INSOLR=0,  IOTSOL=l,  and  ISOLOC=-l.  This  calculation  will  compute  the 
attenuated  solar  plus  thermal  radiation  from  the  top  of  the  atmosphere  to  the  observer  altitude. 
The  result  will  be  written  to  the  file  TAPE24.  The  second  run  of  FASE,  "run  2",  is  set  as  for  the 
downlooking  case:  Hl=Hobs,  H2=Hsurf,  IEMIT=2,  INSOLR=0,  IOTSOL=l,  and  ISOLOC=l. 
This  will  complete  the  calculation  by  using  the  result  from  run  1  as  the  downwelling  radiation 
reaching  Hobs,  which  will  then  be  attenuated  to  Hsurf,  reflected,  and  attenuated  back  to  Hobs. 

4.3.  Line  Coupling 

Line  coupling  is  an  important  phenomenon  that  occurs  in  collisionally  broadened  spectra 
when  the  lines  are  very  close  together.  While  the  theory  will  not  be  discussed  in  this  report,  the 
treatment  of  this  phenomena  in  FASCODE  /  LBLRTM  /  FASE  requires  that  all  of  the  lines  from 
a  particular  line  coupling  region  be  included  in  the  calculation.  The  structure  of  FASCODE  / 
LBLRTM  assumes  that  the  user  recognizes  this  and  adjusts  the  limits  of  the  spectrum 
accordingly.  For  FASE  a  different  strategy  was  adopted  whereby  the  program  checks  to  see  if 
more  lines  need  to  be  included  in  the  calculation  and  adjusts  the  spectral  limits  as  necessary.  To 
avoid  problems  with  output  files,  etc.,  this  adjustment  is  transparent  to  the  user  and  the  output 
file  retains  the  limits  specified  by  the  user.  When  this  occurs  a  warning  message  is  written  to  the 
screen  and  to  the  TAPE6  output  file. 

The  implementation  of  this  limit-checking  scheme  begins  with  the  creation  of  the  TAPE3 
(spectral  line  data  file)  using  the  program  LNFL.  This  program  checks  the  line  coupling 
coefficients  (stored  in  data  statements)  and  puts  the  limits  of  the  various  line-coupling  regions  in 
the  TAPE3  header.  If  the  user  has  selected  a  spectral  range  that  starts  or  ends  in  a  line  coupling 
region,  LNFL  will  modify  the  spectral  limits  to  be  sure  that  all  of  the  coupled  lines  are  selected. 

For  each  atmospheric  layer  FASE  determines  the  spectral  range  for  which  to  retrieve  lines 
from  TAPE3.  After  this  region  is  selected  the  subroutine  CHKLNC  will  check  the  limits  to  see 
if  they  fall  within  a  line-coupling  region.  A  limit  falling  within  a  line-coupling  region  will  be 
reset  to  four  halfwidths  beyond  the  region  (where  the  halfwidth  is  taken  to  be  the  average 
halfwidth  of  the  layer).  Changing  this  limit  will  modify  only  the  spectral  range  of  the  calculation 
and  not  the  spectral  range  of  the  output  file.  This  change  has  been  made  to  both  the  LTE  and 
NLTE  portions  of  the  code  and  the  user  is  notified  in  the  TAPE6  file  if  the  limits  were  changed. 
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4.4.  HITRAN96  Line  Parameters 

A  new  feature  to  FASE  is  that  it  accepts  HITRAN96  line  parameters.  The  code  has  been 
modified  to  accept  36  species  from  the  HITRAN  database  as  opposed  to  32  species,  as  was  the 
case  with  FASCOD3.  The  code  was  modified  by  changing  dimension  statements  and  partition 
function  arrays  to  accommodate  HITRAN96.  It  should  be  noted  that  not  all  of  the  molecules 
have  internally  stored  profiles.  For  species  for  which  there  are  no  internally  stored  profiles  the 
atmosphere  must  explicitly  specified  by  the  user  as  discussed  in  Section  4.5.4. 

4.5.  HITRAN96  Heavy-Molecule  Cross-Sections 

Laboratory  measurements  of  the  absorption  properties  of  heavy  molecules  were  recently 
made  (e.g.  Li  and  Varanasi,  1994;  Varanasi  and  Nemtchinov,  1994,  subsequently  referred  to  as 
"Varanasi's  cross  sections")  and  are  included  on  the  1996  release  of  the  HITRAN  databases 
(Rothman  et  al.,  1992,  1996).  Because  of  the  complexity  of  these  molecules  and  the  fact  that 
individual  spectral  lines  cannot  be  resolved  at  temperatures  and  pressures  typical  of  the  earth’s 
atmosphere,  the  measurements  consist  of  relatively  low  spectral  resolution  information  in  the 
form  of  absorption  cross-sections  (i.e.  without  information  about  line  strengths,  positions, 
halfwidths,  etc.).  Many  of  these  molecules  have  been  included  on  previous  versions  of  the 
HITRAN  databases  in  that  form,  as  temperature  dependent  cross  sections  rather  than  the 
traditional  HITRAN  format  of  line  position,  line  strength,  halfwidth,  etc.  However,  the  Varanasi 
measurements  are  provided  as  a  function  of  pressure  and  temperature  as  opposed  to  the  previous 
HITRAN  format,  which  extrapolated  the  temperature-dependent  values  to  zero  pressure.  To 
compute  the  optical  depth  for  a  particular  (p,T)  combination  with  the  'old'  data,  it  was  necessary 
to  interpolate  between  temperatures  while  convolving  with  the  appropriate  Lorentz  line  shape  for 
the  pressure,  finally  combining  this  derived  cross  section  with  the  density  weighted  path 
amounts.  With  the  Varanasi  data,  one  is  required  to  interpolate  and/or  extrapolate  between  the 
tabulated  data  to  arrive  at  the  correct  absorption  coefficient.  Thus,  in  order  to  accommodate  this 
new  data  directly,  it  would  be  necessary  to  re-configure  the  FASE  algorithms  to  bypass  the 
convolution  and  perform  the  necessary  interpolation  and/or  extrapolation.  Further,  one  must 
develop  an  appropriate  scheme  for  extrapolating  the  tabulated  values  to  pressures  and 
temperatures  outside  the  range  of  the  measurements. 

One  approach  through  which  the  Varanasi  cross-sections  may  be  utilized  by  the  FASE 
algorithms  is  to  perform  a  least-squares  fit  to  the  data  with  a  series  of  spectral  lines,  solving  for 
the  line  strength,  halfwidth,  and  lower  state  energy  (Toon,  G.,  1995).  The  pseudo-line  approach 
incorporates  all  the  (p,T)  combinations  of  data  into  the  least-squares  fit  and  increases  the 
accuracy  from  a  simple  interpolation  between  two  (p,T)  values.  The  pseudo-line  data  is  easily 
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merged  with  the  spectral  data  taken  from  HITRAN  when  creating  the  TAPE3’  line  file,  and  the 
radiative  transfer  calculations  are  performed  as  usual,  with  the  model  using  the  spectral  line 
routines,  rather  than  the  cross-section  routines,  when  computing  the  absorption  of  these  species. 

One  drawback  to  the  pseudo-line  approach  is  that  new  sets  of  pseudo-lines  are  required  as 
more  molecules  are  added.  In  the  long  run  it  would  be  desirable  to  develop  a  scheme  to  use  the 
cross-sections  given  as  a  function  of  pressure  and  temperature.  Since  the  problem  lies  in 
interpolating/extrapolating  the  absorption  coefficient  data  to  specific  values  of  the  pressure  and 
temperature,  one  might  use  the  pseudo-line  approach  to  fill  out  the  full  matrix  of  required 
pressures  and  temperatures  (that  is,  use  FASE  to  compute  the  absorption  coefficients  for  the 
regions  not  represented  by  the  measurements).  For  the  time  being  it  is  sufficient  to  investigate 
the  implementation  of  the  pseudo-line  approach  and  begin  to  explore  ways  in  which  the  actual 
data  can  be  used,  perhaps  in  conjunction  with  'pseudo-data'  based  on  calculations  with  the 
pseudo-lines. 

The  pseudo-line  approach  has  been  evaluated  using  line  data  for  CFC-12  from  fits  made  by 
Toon  (1995).  Before  this  method  can  be  accepted  there  are  two  questions  that  must  be 
addressed:  (1)  does  the  calculation  with  the  pseudo-lines  accurately  reproduce  the  measured 
absorption  coefficients?  (2)  what  is  the  timing  impact  of  using  the  pseudo-lines  compared  to  the 
(low  spectral  resolution)  cross-sections?  These  are  discussed  in  the  following  two  sections. 
Section  4.5.3  outlines  the  coding  changes  required  to  implement  the  pseudo-line  spectral 
information  while  Section  4.5.4  explains  how  one  can  do  a  calculation  utilizing  this  data. 

4.5.1 .  Accuracy  of  Pseudo-Line  Approach 

The  accuracy  of  the  pseudo-line  approach  was  evaluated  by  computing  the  optical  depth  for  a 
horizontal  path  with  constant  pressure  and  temperature  and  a  single  molecular  species  (CFC-12). 
The  absorption  coefficient,  k,  can  be  determined  directly  from  the  calculation: 


where  x  is  the  optical  thickness,  N  is  the  absorber  density  integrated  along  the  line-of-sight  path 
and  k  has  units  of  cm2.  The  calculated  absorption  coefficient  can  then  be  compared  directly  with 
the  measurements  in  order  to  determine  the  validity  of  this  approach. 

This  test  was  done  for  several  different  combinations  of  pressure  and  temperature  as  well  as 
different  absorption  path  lengths.  Figure  7  and  Figure  8  show  the  calculated  absorption  cross- 
section  for  F-12  for  two  different  (p,T)  combinations  and  a  1  km  path.  The  percent  difference 
between  the  computed  cross  sections  and  Varanasi's  measurements  are  plotted  in  Figure  9 
through  Figure  12  for  the  1  km  and  50  km  paths.  Because  of  systematic  differences  in  the  wings 
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of  the  band,  where  the  absorption  is  very  small,  only  the  region  from  860  -  940  cm'1  is  shown  in 
the  difference  plots  (note  that  the  scales  are  different  between  Figure  9/Figure  10  and  Figure 
1 1/Figure  12).  The  systematic  differences  are  clearly  seen  in  Figure  13,  a  plot  of  the  difference 
rather  than  percent  difference,  and  Figure  24,  where  it  can  be  seen  that  negative  values  for  the 
data  were  set  to  zero.  It  is  interesting  to  note  that  the  percent  difference  over  the  region  plotted  in 
Figure  9  through  Figure  12  is  identical  regardless  of  the  length  of  the  absorption  path.  The  fact 
that  the  percent  difference  does  not  scale  with  path  length  indicates  a  consistency  in  the  FASE 
calculation.  Figure  14  through  Figure  18  detail  certain  regions  of  Figure  7,  while  Figure  19 
through  Figure  23  detail  regions  of  Figure  8.  In  these  figures  the  computed  value  of  the  cross 
section  (solid  line)  is  compared  with  the  data  measured  and  reported  by  Varanasi  (dotted  line). 
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Figure  7  Absorption  Cross-Section  from  the  FASE  Calculation  for  a  1km  Path  at  p=169.9 
torr,  T=216.0  K  (path  amount  of  F-12=1.814E+14  molecules  cm-2) 
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Figure  8  Absorption  Cross-Section  as  Computed  from  the  FASE  Calculation  for  a  1km 
Path  at  p=700.3  torr,  T=296.2  K  (path  amount  of  F-12=5.437E+14  molecules  cm"2) 


Figure  9  Percent  Difference  Between  the  Cross-Section  Computed  with  Pseudo-Lines  and 
the  Data  from  Varanasi  for  a  1km  Path  at  p=169.9  torr,  T=216.0  K 
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p  =  1  69.9  torr,  T=216.0  K,  50km  path 


wavenumber  (cm"’) 


Figure  10  Percent  Difference  Between  the  Cross-Section  Computed  with  Pseudo-Lines  and 
the  Data  from  Varanasi  for  a  50km  Path  at  p=169.9  torr,  T=216.0  K 


p  =  700.3  torr,  1=296.2  K,  1km  path 
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Figure  11  Percent  Difference  Between  the  Cross-Section  Computed  with  Pseudo-Lines  and 
the  Data  from  Varanasi  for  a  1km  Path  at  p=700.3  torr,  T=296.2  K 
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Figure  12  Percent  Difference  Between  the  Cross-Section  Computed  with  Pseudo-Lines  and 
the  Data  from  Varanasi  for  a  50km  Path  at  p=700.3  torr,  T=296.2  K 


p=169.9  torr,  T=216  K,  1km  path,  solid  =  FASE 
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Figure  13  Difference  Between  Computed  and  Measured  Cross-Section  for  a  1km  Path  at 

p=169.9  torr,  T=216.0  K 
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p=169.9  to rr,  T=216  K,  1km  path,  solid  —  FA5E 


Figure  14  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  7  Over  the  Range  860-880  cm 1 
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Figure  15  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  7  Over  the  Range  880-900  cm 1 
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Figure  16  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  7  Over  the  Range  900-910  cm'1 


Figure  17  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  7  Over  the  Range  920-925  cm'1 
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p-169.9  to rr,  T=216  K,  1km  path,  solid  =  FASE 
0  050  f=  :  T  1  T  1  T - T - T - - - 1 - T - 11 - T - T - i - 1 - ' — 


0.040  r 


o  p 

7  0.030  f- 

'E  E 


D 


0.01  0  L _ * _ 1  ■  I  :  .  1  i  1  I  I  !  ■  I  I  I  I  ■  .4 

942.0  9^2.5  943.0  943.5  944.0 

wcvenumber  (cm-1) 


Figure  18  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  7  Over  the  Range  942-944  cm'1 

p=700.3  tor-,  T=296.2K,  1km  poth,  solid  =  FASE 
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Figure  19  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  8  Over  the  Range  860-880  cm'1 
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Figure  20  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  8  Over  the  Range  880-900  cm'1 


Figure  21  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  8  Over  the  Range  900-910  cm 1 
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p  =  700.3  torr,  7=296. 2K,  1km  path,  solid  =  FASE 


Figure  22  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured  Cross- 
Section  (Dotted  Line)  for  the  Conditions  Given  in  Figure  8  Over  the  Range  920-925  cm 1 
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Figure  23  Comparison  of  the  Computed  Cross-Section  (Solid  Line)  with  the  Measured 
Cross-Section  (Dotted)  for  the  Conditions  Given  in  Figure  8  Over  the  Range  942-944  cm 1 
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CCI2F2  T,p:  216.0000  169.9 


Figure  24  Measurements  of  F-12  in  the  Far  Wings  of  the  Band 

In  general,  there  is  a  slight  offset  between  the  calculated  cross-sections  and  the  Varanasi  data 
for  the  169.9  torr,  216.0  K  case  (Figure  14  through  Figure  18).  The  peak-to-peak  magnitude  of 
the  features  in  Figure  14  is  somewhat  less  for  the  calculated  cross-sections  than  for  the  data.  In 
Figure  15  the  peak-to-peak  values  are  better,  but  there  is  some  difference  in  the  shape  of  the 
calculation.  This  is  evident  around  887,  893,  and  895-896  cm1.  The  differences  in  shape  are 
especially  noticeable  in  Figure  16,  from  900  -  910  cm'1,  which  has  the  worst  agreement  for  this 
case  (except  for  the  wing  regions  where  the  absorption  is  very  weak).  The  region  around  the 
band  center  is  plotted  in  Figure  17  and  one  can  see  that  there  is  a  slight  offset,  but  otherwise  very 
good  agreement  between  the  computation  and  the  data.  Finally,  Figure  18  shows  a  region  in  the 
wings  of  the  band  where  there  appears  to  be  much  more  noise  than  real  data. 

The  comparisons  of  the  measured  data  with  calculations  using  the  pseudo-lines  for  770.3  torr 
and  296.2  K  are  given  in  Figure  19  through  Figure  23.  These  figures  show  a  much  more  distinct 
offset  than  Figure  14  through  Figure  18.  While  the  structure  is  difficult  to  determine  in  Figure 
18,  it  appears  to  match  the  data  very  well  (especially  from  877  -  878  cm1).  The  peak-to-peak 
agreement  is  also  very  good  in  Figure  19  and  Figure  21,  while  Figure  20  clearly  shows  the  offset. 
As  with  Figure  18,  Figure  23  gives  an  indication  that  the  agreement  in  the  wings  of  the  band  is 
very  poor. 
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The  above  comparisons  lead  to  another  aspect  of  the  problem  that  must  be  evaluated,  namely 
the  determination  of  meaningful  values  for  the  cross-section  and  the  establishment  of  a  frequency 
cutoff  for  the  inclusion  of  pseudo-lines  in  the  calculation.  What  is  necessary  is  a  reasonable 
criterion  for  determining  a  wavenumber  cutoff  whereby  the  noise  in  the  data  coupled  with 
uncertainties  in  the  calculation  combine  to  make  it  impossible  to  determine  the  validity  of  the 
cross-section.  What  will  be  important  is  a  combination  of  the  small  magnitude  of  the  cross- 
section  in  this  cutoff  region  and  the  magnitude  of  overlapping  contributions  of  other  species 
(which  may  completely  mask  the  weak  absorption  of  the  cross-section).  A  means  of 
understanding  the  errors  in  the  cross  section  computed  using  the  pseudo-lines  and  in  the  Varanasi 
data  can  be  obtained  through  a  propagation  of  errors  analysis. 

As  described  above,  the  absorption  coefficient  computed  by  FASE  is  the  optical  thickness 
divided  by  the  path  absorber  amount.  Thus  the  error  in  the  absorption  coefficient  can  be 
expressed  as  the  sum  of  the  individual  errors: 


Akcalc  =  All 
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Eq.  20 


which  reduces  to 


Ak„lt=2(AT-ANkciJ 


Eq.  21 


For  the  calculations  shown  above,  the  appropriate  values  for  Eq.  21  are  given  in  Table  14.  These 
values  yield  an  uncertainty  in  the  absorption  coefficient  of  about  1%  for  k=1.0e-17  (center  of 
band),  and  an  uncertainty  approaching  10%  in  the  wings  of  the  band. 

The  errors  for  the  measured  data  can  be  determined  in  two  ways.  First,  a  qualitative 
examination  of  the  measurements  in  the  wings  of  the  band  shows  a  random  variation  of  about 
0.1%.  A  somewhat  more  quantitative  method  is  to  use  the  reported  error  bars  (Varanasi  and 
Nemtchinov,  1994)  for  the  variation  in  the  integrated  cross-section:  AS=0.07e-17,  where  S  is  the 
integral  of  k  over  the  entire  band: 


O  f  k™as(c)da  _  Y  kmeas(c>)AG 

J  N  ^  N 

o  o 


Eq.  22 


and  N0  is  used  to  convert  the  units  of  k  (cm'1  atm'1)  to  the  given  units  of  S  (cm  molecule'1). 

Note  that  the  units  are  slightly  different  from  the  ones  used  in  earlier  equations,  but  they  conform 
to  the  reported  values.  If  one  assumes  that  AS  and  N0  are  known,  the  variation  in  S  can  be 
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attributed  to  the  variations  in  k.  Assuming  that  the  variations  in  k  are  random  throughout  the 
band,  the  error  in  k  can  be  written  as: 


Akmeas(a)  = 


AS(q)No 

5>a 


Eq.  23 


where  the  denominator  represents  the  entire  bandwidth.  Using  the  values  in  Table  14,  the  error 
in  the  absorption  coefficient  is  less  than  0.2%,  which  is  similar  to  that  determined  by 
examination  of  the  measurement.  In  this  table,  N  is  the  path  absorber  amount  (molecules  cm'2) 
and  the  1%  error  is  attributed  to  how  this  number  is  written  to  the  output  files;  x  is  the  optical 
thickness  and  its  error  is  the  stated  numeric  error  of  FASE;  k  is  the  absorption  cross-section 
(cm2),  AS  is  the  error  in  the  integrated  cross-section  (cm  molecule*1)  (from  Varanasi  and 
Nemtchinov,  1994);  N0  is  Loschmidt's  number  (molecule  cm-3  atm*1);  and  Act  is  the  spectral 

band  width  (cm*1).  Because  of  the  weak  signal  in  the  wings  of  the  band,  and  because  the 
assumption  that  the  error  in  the  integrated  cross-section  can  be  attributed  to  the  absorption 
coefficient  equally  throughout  the  band,  it  is  not  clear  if  either  of  these  error  values  are 
appropriate  when  examined  individually.  However,  their  consistency  would  indicate  that  0.1- 
0.2%  represents  a  reasonable  guess. 

Table  14  Values  Used  in  Propagation  of  Error  Calculations 


Quantity 

Value 

(order  of  magnitude) 

N 

1.0e+14 

AN 

1%  =  1.0e+12 

X 

1.0e-03 

Ax 

0.5%  =  5.0e-06 

k 

1.0e-17 

AS 

0.013e-17 

N0 

2.687e+19 

A  a 

100 

4.5.2.  Pseudo-Line  Tinning  Tests 

The  ability  to  read  the  new  cross-section  data  as  a  series  of  pseudo-lines  is  an  inexpensive 
alternative  to  developing  a  more  sophisticated  interpolation/extrapolation  scheme.  However,  the 
price  paid  is  one  of  algorithm  timing.  The  inclusion  of  pseudo-lines  can  add  anywhere  from 
thousands  to  tens  of  thousands  of  spectral  lines  to  a  given  calculation.  For  calculations,  which 
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require  a  large  number  of  atmospheric  layers,  this  can  add  a  significant  amount  of  time  to  the 
calculation. 

A  series  of  tests  was  conducted  in  the  spectral  range  from  725  -  1245  cm'1  using  the  heavy 
molecule  species  F-l  1,  F-12,  and  CC14.  These  species  were  chosen  because  they  have  data  in  the 
form  of  'old'  cross-sections  and  'new'  pseudo-lines.  Calculations  were  done  for  the  uplooking 
radiance  and  transmittance  from  0  -  40  km  and  0  -  100  km  with  a  US  Standard  Atmosphere.  The 
timing  results  given  in  Table  15  represent  an  average  of  eight  individual  runs  for  each  case,  and 
the  standard  deviation  for  each  case  was  less  than  3  seconds. 

The  calculations  from  100  km  to  the  surface  take  longer  than  those  from  40  km  because  of  an 
increased  number  of  atmospheric  layers,  and  because  the  upper  layers  require  higher  spectral 
resolution  than  the  lower  layers.  With  the  cross-sections,  calculations  from  100  km  take 
approximately  3.1  times  longer  than  those  from  40  km.  The  difference  is  larger  with  the  pseudo¬ 
lines,  a  factor  of  about  3.5,  because  each  pseudo-line  must  be  accurately  computed  at  each  layer. 
In  reality,  the  concentration  profiles  of  many  species,  the  heavy  molecules  in  particular,  decrease 
rapidly  with  altitude.  Using  the  option  to  zero  the  contribution  from  weak  lines  when  they  no 
longer  contribute  significantly  to  the  overall  optical  depth  (NOZERO=0  on  input  RECORD  3.1) 
dramatically  decreases  the  computation  time  for  both  the  pseudo-line  and  cross-section  cases. 
For  the  0-100  km  case  shown  here  there  is  a  20%  timing  penalty  for  using  the  pseudo-lines 
instead  of  the  cross-sections,  which  is  small  compared  to  the  change  in  accuracy  with  this  new 
data. 


Table  15  Timing  Test  Results 


Case 

Old  Cross-Section 

New  Pseudo-Line 

0-40 

178.96 

401.11 

0-  100 

562.93  J 

1413.07 

0-100 

(zero  small 

absorber  amounts) 

263.02 

327.36 

4.5.3.  Pseudo-Line  Algorithm  Implementation 
4.5. 3.1.  Line  File  Creation  (LNFL) 

A  minimal  amount  of  changes  were  required  to  LNFL  (the  program  used  to  select  line 
parameters  from  the  HITRAN  database  and  put  them  in  the  correct  format  for  FASCODE  and 
LBLRTM)  in  order  to  accommodate  the  pseudo-lines.  The  lines  are  read  from  an  'external  data 
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tape'  and  merged  with  the  lines  from  HITRAN.  There  are  36  molecules  in  the  current  version  of 
HITRAN  and  the  pseudo-lines  were  given  IDs  starting  with  51.  Consequently,  the  dimensions  of 
arrays,  the  parameters  "NTMOL"  and  "NSPECI",  and  read/write  format  statements  needed  to  be 
changed  to  accommodate  the  additional  molecules.  Data  statements  relating  to  the  molecule 
name,  the  number  of  isotopes  (1)  and  the  isotope  ED  number  (assumed  to  be  111)  were  also 
changed,  and  the  temperature  dependence  of  the  halfwidth  was  set  at  0.5  (Toon,  1995). 

4.5.3.2.  Main  Program  and  Subroutines 

In  all  of  the  files  listed  below,  the  parameter  statements  "NTMOL"  (the  total  number  of 
molecules  allowed)  and  "NSPECI"  (the  total  number  of  isotopes  allowed)  were  changed  from  36 
and  79  to  64  and  90.  These  parameters  are  used  to  define  array  sizes  in  dimension  statements; 
additional  species  will  require  further  changes  to  "NSPECI",  but  not  "NTMOL"  (until  the  total 
number  of  molecules  is  greater  than  50). 


fasatm.f 

The  name  of  the  molecules  ('HMOLEC')  and  the  corresponding  molecular  weights 
("AERMWT")  were  set  appropriately.  The  molecular  weight  was  set  to  1.0  for  each  of  the 
pseudo-line  species  so  that  the  Doppler  width  will  be  essentially  the  same  as  the  grid  spacing  and 
the  pseudo-lines  can  never  become  too  strong  (Toon,  1995). 

faseOl.f  /  nonlte.f 

The  only  changes  are  to  the  parameter  statements  listed  above. 

oprop.f 

In  addition  to  the  parameter  statements  listed  above  and  molecular  weight  ("SMASSI")  and 
isotope  information  ("ISONM"  and  "IS082"),  the  changes  to  this  file  involve  the  partition 
functions  and  degeneracy  factors.  The  inputs  required  for  calculating  the  partition  sums  were 
provided  by  Toon  (1995).  The  formulas  for  the  vibrational  partition  function  (Qv),  the 
stimulated  emission  term  (Sc),  and  the  rotational  partition  function  (Qr)  are  given  as: 
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Eq.  24 


with  b  -  1.5,  the  fundamental  vibrational  frequencies,  Dj,  listed  in  Table  16,  and  1),  the  center 
frequency  of  the  line  in  question. 

Table  16  Fundamental  Vibrational  Frequencies  for  F-12 


vi 

(cm1) 

1102 

667 

468 

262 

321 

1161 

462 

923 

437 


4.5.4.  Execution  of  FASE  with  Pseudo-Lines 

Atmospheric  input  to  FASE  is  accomplished  through  the  use  of  the  module  "FASATM". 
"F ASATM"  is  used  to  select  the  atmospheric  profile  to  be  used  by  the  module  and  sets  up  the 
layers.  "FASATM"  is  set  from  the  "IATM"  parameter,  contained  in  RECORD  1.2.  When 
IATM  is  set  to  1,  FASE  is  instructed  to  look  to  RECORD  3  for  specification  of  the 
atmospheric  profile.  Input  parameters  for  "FASATM”  are  located  in  RECORDS  3.1  through  3.6 
of  the  TAPES  input  file. 

The  number  of  molecular  species  desired  is  set  from  RECORD  3.1  with  the  "NMOL" 
parameter.  NMOL  must  be  set  to  be  at  least  equal  to  the  highest  number  of  the  molecule 
desired.  Molecules  from  the  HITRAN  Database  have  been  given  molecule  identification 
numbers  1-36  while  pseudo-lines  are  given  molecule  identification  numbers  51-60.  Some  of  the 
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HITRAN  96  and  pseudo-line  molecules  do  not  have  internally  stored  profiles  (see  Table  17)  and 
for  these  species  the  user  must  explicitly  specify  the  atmosphere. 

FASE  employs  a  user-defined  atmospheric  profile  when  the  parameter  "MODEL"  is  set  to  0 
(zero)  in  RECORD  3.1.  RECORDS  3.4  through  3.6  contain  the  information  required  for  a  user- 
defined  atmospheric  profile.  Contained  within  RECORD  3.4,  "IMMAX"  is  the  number  of 
atmospheric  profile  boundaries  to  be  read  in  by  the  model.  RECORD  3.5  contains  the  physical 
parameters  required  to  define  a  boundary. 

Table  17  Molecules  Without  Profiles 


HITRAN-96 

Pseudo-Lines 

ID  Number 

Molecule 

ID  Number 

Molecule 

29 

COF, 

52 

CF„ 

30 

SF, 

57 

F-142b 

31 

H,S 

60 

SF, 

32 

HCOOH 

33 

HO, 

34 

O-Atom 

Note  that  if  one  wishes  to  run  with  the  first  seven  molecules  and  pseudo-lines,  one  can 
specify  a  profile  for  molecules  1-7  and  the  pseudo-line,  and  zero  for  the  remainder  of  the  species. 
Also  note  that  FASE  does  not  check  to  see  if  more  than  one  of  the  same  species  has  been 
selected.  For  example,  CL0N02  has  lines,  pseudo-lines,  and  cross-sections. 

4.6.  Cloud/Rain  Upgrades 

Changes  to  the  cloud  and  rain  routines  have  been  employed  based  on  improvements  that  have 
been  developed  for  MODTRAN  (Berk,  1995).  This  set  of  upgrades  makes  it  easier  to  modify  the 
cloud/rain  models  (by  re-configuring  and  clarifying  the  coding  structure)  as  well  as  providing  for 
a  more  flexible  prescription  of  realistic  cloud  and  aerosol  layers  and  their  associated  optical 
properties.  The  improved  layering  scheme  allows  for  multiple  overlapping  and  non-overlapping 
clouds.  For  cumulus  and  stratus  type  clouds,  with  and  without  rain,  the  upgrades  include: 

•  adjustable  cloud  parameters  (thickness,  altitude,  vertical  extinction,  water/ice  column 
amounts,  humidity,  and  scattering  phase  functions), 

•  decoupling  of  the  clouds  from  aerosols  (allowing  clouds  and  aerosols  to  be  present  at  the 
same  altitude), 
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•  the  introduction  of  ice  particles,  and 

•  a  flexible  set  of  user-defined  cloud  spectral  properties,  and  water  droplet,  ice  particle,  and 
rain  rate  profiles. 

Cloud  profiles  are  merged  with  other  atmospheric  profiles  by  combining  and  adding  layer 
boundaries  where  necessary.  The  cloud  profiles  will  be  merged  with  existing  layers  if  the  cloud 
layer  boundary  is  within  0.5m  of  an  atmospheric  boundary  level,  otherwise  a  new  atmospheric 
boundary  level  is  created.  (The  number  of  available  levels  is  controlled  by  the  parameter 
"LAYDIM").  The  cloud/rain  profiles  used  by  the  code  are  written  to  TAPE6. 

Figure  25  illustrates  calculations  with  the  new  cloud  models.  The  cloud  consisted  of  a 
cumulus  cloud  with  overlying  cirrus  (the  cloud  top  and  cirrus  thickness  were  held  constant)  and 
the  figure  shows  the  change  in  radiance  with  a  change  in  cloud  thickness  and  base  height.  The 
marked  change  in  radiance  for  a  cloud  with  a  base  at  0.4  km  and  one  at  1.88  km  shows  the  need 
for  a  versatile  cloud  simulation  capability,  even  for  direct  radiance  and  transmittance  (no 
scattering)  calculations. 


Varying  Cioud  Base  and  thickness 


Figure  25  Example  of  the  New  FASE  Cloud  Options:  User  has  the  Ability  to  Change 

Cloud  Base  Altitude  and  Thickness 
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5.  Suggestions  for  Future  Upgrades  and  ModelVaiidation 

There  are  several  areas  of  development  in  which  FASE  is  still  incomplete.  Of  these,  the  most 
significant  is  the  absence  of  a  capability  for  including  multiple  scattering.  In  addition,  there  are  a 
number  of  validations  that  should  be  performed  in  order  to  confirm  the  continued  accuracy  of  the 
model.  The  following  section  discusses  suggested  areas  for  future  development. 

5.1.  Multiple  Scattering 

Multiple  scattering  was  eliminated  from  the  development  of  LBLRTM.  For  FASE  it  was 
decided  that  it  was  not  cost-effective  to  re-implement  the  2-stream  multiple  scattering  routines 
that  had  been  a  part  of  FASCODE.  While  FASE  does  have  the  option  to  create  the  output  files 
that  can  be  used  as  input  to  multiple  scattering  codes,  a  coherent  interface  between  FASE  and 
one  such  code  should  be  pursued.  The  likely  candidate  for  coupling  with  FASE  is  the  multiple 
scattering  model  CHARTS.  The  advantage  to  this  code  (over  the  other  major  code,  DISORT)  is 
that  it  is  designed  to  work  with  the  panel  structure  of  LBLRTM,  rather  than  the  monochromatic 
input  required  for  DISORT.  That  is,  the  amount  of  calculation  overhead  is  reduced  when 
running  CHARTS  over  a  broad  spectral  region,  making  it  uniquely  suited  for  coupling  with 
FASE. 

5.2.  Cloud/Rain  Upgrades 

There  are  now  further  upgrades  to  the  cloud  and  rain  routines  available  from  MODTRAN4. 
To  maintain  consistency  in  the  common  elements  of  the  Air  Force  radiance  and  transmittance 
codes,  these  upgrades  should  be  added  to  FASE. 

5.3.  Beta-Test  Results 

FASE  has  not  been  released  to  the  general  public,  but  rather  to  a  limited  set  of  users  (‘beta- 
testers  )  who  were  encouraged  to  report  problems  with  the  code  and  suggest  new  features  or 
enhancements  to  the  existing  code.  While  most  of  the  problems  have  been  fixed,  a  few  issues 
remain,  mainly  related  to  the  new  cloud  and  aerosol  routines.  Beta-test  users  of  FASE  at  the  Air 
Force  Research  Laboratory-Hanscom  have  identified  problems  in  the  new  cloud/aerosol  routines 
when  doing  calculations  in  the  microwave.  Investigation  of  the  cause  of  these  problems  led  to 
the  identification  of  several  errors  in  the  implementation  of  the  upgrades  from  MODTRAN. 
Fixing  these  errors  did  not  solve  the  microwave  problems,  and  this  work  should  continue. 


56 


5.4.  Calibration  and  Validation 

An  important  part  of  the  acceptance  of  FASCODE/FASE  over  other  publicly  available 
radiative  transfer  models  is  the  on-going  calibration  and  validation.  Examples  of  this  include 
participation  in  the  ICRCCM  testing  and  the  validation  of  the  model  with  ARM  data  . 
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6.  Summary 


Significant  progress  has  been  made  in  the  merging  of  FASCODE  and  LBLRTM  to  create 
FASE.  The  result  is  a  radiative  transfer  model  that  contains  state-of-the-art  atmospheric  physics 
through  validations  of  LBLRTM  with  ARM  data,  while  also  incorporating  many  of  the  features 
from  FASCODE  which  are  not  found  in  LBLRTM  but  are  required  for  flexible  use  of  the  codes 
under  airborne,  ground-  and  space-based  conditions.  Future  updates  to  FASE  should  include  a 
coupling  of  the  code  to  multiple-scattering  routines.  Further,  FASE  should  continue  to  undergo 
comparisons  with  other  models  and  validation  against  measurements  from  the  DOE  ARM  site  as 
well  as  a  variety  of  DoD  datasets. 


58 


7.  References 


Acharya,  P.K.,  D.C.  Robertson,  and  A.  Berk,  Upgraded  Line-of-Sight  Geometry  Package  and 
Band  Model  Parameters  for  MODTRAN,  PL-TR-93-2127,  Phillips  Laboratory 
Directorate  of  Geophysics,  Hanscom  AFB,  Massachusetts,  May  1993.  ADA280952. 

Anderson,  G.  P.,  F.X.  Kneizys,  E.P.  Shettle,  L.W.  Abreu,  J.H.  Chetwynd,  R.  E.  Huffman,  and 
L.  A.  Hall,  UV  Spectral  Simulations  Using  LOWTRAN7,  in  Advisory  Group  for 
Aerospace  Research  and  Development  (AGARD)  Conference  Proc.  454,  1990. 

Anderson,  G.  P.,  et  al„  History  of  One  Family  of  Atmospheric  Radiative  Transfer  Codes,  The 
European  Symposium  on  Satellite  Remote  Sensing,  Conference  on  Passive  Infrared 
Remote  Sensing  of  Clouds  and  Atmosphere  II,  26-30  September  1994,  Rome,  Italy. 

Anderson,  G.  P.„  et  al„  1995,  FASCODE/MODTRAN/LOWTRAN:  Past/Present  Future, 
Proceedings  of  the  18lh  Annual  Review  Conference  on  Atmospheric  Transmission 
Models,  6-8  June,  1995,  Phillips  Lab/Geophysics  Directorate,  Hanscom  AFB. 

Anderson,  G.  P.,  Phillips  Laboratory  Directorate  of  Geophysics,  Private  Communication, 
1994. 

Berk,  A.,  Upgrades  to  the  MODTRAN  Layer  Cloud/Rain  Models,  Rpt.  No.  SSI-SR-56, 
Spectral  Sciences,  Inc.,  99  S.  Bedford  St.,  Burlington,  MA,  01803,  1995. 

Chetwynd,  J.H.,  Phillips  Laboratory  Directorate  of  Geophysics,  Private  Communication, 
1994. 

Clough,  S.A.  and  F.X.  Kneizys,  Convolution  Algorithm  for  the  Lorentz  Function,  Applied 
Optics,  18,  2329,  1979. 

Clough,  S.A.,  F.X.  Kneizys,  L.  S.  Rothman,  and  W.O.  Gallery,  Atmospheric  spectral 
transmittance  and  radiance:  FASCODE1B,  SPIE  V.277,  Atmospheric  Transmission 
1981. 

Clough,  S.A.,  et  al..  Radiative  Transfer  Model  Development  in  Support  of  the  Atmospheric 
Radiation  Measurement  (ARM)  Program,  Proceedings  of  the  Second  ARM  Science  Team 
Meeting,  DOE  Conf.  9303112,  Norman,  OK,  1-4  March  1993. 

Clough,  S.A.  and  P.  Brown,  Atmospheric  and  Environmental  Research,  Inc,  Private 
Communication,  1994. 

Edlen,  K.,  The  Refractive  Index  of  Air,  Metrologia,  2,  12,  1966. 

Gallery,  W.O.,  F.X.  Kneizys,  and  S.A.  Clough,  Air  Mass  Computer  Program  for 
Atmospheric  Transmittance/Radiance  Calculations:  FSCATM,  AFGL-TR-83-0065,  Air 
Force  Geophysics  Laboratory,  Hanscom  AFB,  Massachusetts,  (NTIS  AD  A  132108) 
1983. 


59 


Gallery,  W.O.  and  S.A.  Clough,  FFTSCAN:  A  Program  for  Spectral  Smoothing  Using 
Fourier  Transforms,  PL-TR-92-2131,  Geophysics  Directorate/Phillips  Lab,  Hanscom 
AFB,  Massachusetts,  1992. 

Hall,  L.A  .,  and  G.  P.  Anderson,  High-Resolution  Solar  Spectrum  between  2000  and  3100  A, 
J.  Geophys.  Res.,  96,  12927-12931,  1991 . 

Isaacs,  R.G.  W.-C.  Wang,  R.  D.  Worsham  and  S.  Goldenberg,  Multiple  Scattering 
LOWTRAN  and  FASCODE  Models,  Applied  Optics,  26, 1272-1281, 1987. 

Li,  Z.,  and  P.  Varanasi,  Measurement  of  the  Absorption  Cross-Section  of  CFC-11  at 
Conditions  Representing  Various  Model  Atmospheres,  J.  Quant.  Spectrosc.  Radiat. 
Transfer,  52,  137-144, 1994. 

Kneizys,  F.X.,  E.P.  Shettle,  W.O.  Gallery,  J.H.  Chetwynd,  Jr.,  L.W.  Abreu,  J.E.A.  Selby, 

R. W.  Fenn,  and  R.A.  McClatchey,  Atmospheric  Transmittance/Radiance:  Computer 
Code  LOWTRAN  5,  AFGL-TR-80-0067,  Air  Force  Geophysics  Laboratory,  Hanscom 
AFB,  Massachusetts,  21  February  1980.  ADA088215. 

Kneizys,  F.X.,  E.P.  Shettle,  W.O.  Gallery,  J.H.  Chetwynd,  Jr.,  L.W.  Abreu,  J.E.A.  Selby, 

S. A.  Clough,  and  R.W.  Fenn,  Atmospheric  Transmittance/Radiance:  Computer  Code 
LOWTRAN  6,  AFGL-TR-83-0187,  Air  Force  Geophysics  Laboratory,  Hanscom  AFB, 
Massachusetts,  1  August  1983.  ADA137689. 

Kurucz,  R.L.,  The  Solar  Irradiance  by  Computation,  in  Proc.  of  the  17th  Annual  Review 
Conference  on  Atmospheric  Transmission  Models,  7  June  1994,  Geophysics 
Directorate/Phillips  Laboratory,  Hanscom  AFB,  MA,  1994. 

Miller,  S.M.,  H.  E.  Snell,  and  J.-L.  Moncet,  Simultaneous  Retrieval  of  Middle  Atmospheric 
Temperature  and  Trace  Gas  Species  Volume  Mixing  Ratios  from  CIRRIS,  J.  Geophys. 
Res.,  104, 18,697-18,714,  1999. 

Minschwaner,  K.,  G.  P.  Anderson,  L.  A.  Hall,  and  K.  Yoshino,  Polynomial  Coefficients  for 
Calculating  O2  Schumann-Runge  Cross  Sections  at  0.5  cm-1  Resolution,  J.  Geophys. 
Res.,  97,  10103-10108, 20  June  1992. 

Minschwaner,  K.,  G.  P.  Anderson,  L.  A.  Hall,  R.  J.  Thomas,  D.  Rusch,  A.  Berk,  and  J. 
Conant,  Scattered  Ultraviolet  Radiation  in  the  Upper  Stratosphere  3:  Modeling  and 
Analysis,  J.  Geophys.  Res.,  in  press  1995. 

Moncet,  J.-L.,  and  S.A.  Clough,  CHARTS:  Code  for  high  resolution  accelerated  radiative 
transfer  with  scattering,  in  IRS  '92:  Current  Problems  in  Atmospheric  Radiation, 
Keevallik  and  Kamer,  editors,  A.  Deepak,  pp.  493-494,  1992. 

Ridgway,  W.L.,  R.A.,  Moose,  and  A.  C.  Cogley,  Atmospheric  Transmittance/Radiance 
Computer  Code  FASCOD2,  Air  Force  Geophysics  Laboratory  Technical  Report  AFGL- 
TR-82-0392,  1982.  ADA130241. 

Rothman,  L.  S.,  et  al..  The  HITRAN  Molecular  Database:  Editions  of  1991  and  1992, 
J.Q.S.R.T.,  48, 469-507,  1992. 


60 


Rothman,  L.  S.,  and  A.  McCann,  HITRAN96,  Phillips  Laboratory  Geophysics  Directorate, 
Hanscom  AFB,  Massachusetts,  released  on  CD-ROM,  1996. 

Shettle,  E.P.,  and  S.M.  Anderson,  Ozone  Cross-Sections  for  the  Chappuis  and  Wulf 
Absorption  Bands,  private  communication  from  G.P.  Anderson,  1995. 

Smith,  H.J.P.,  D.  J.  Dube,  M.E.  Gardner,  S.A.  Clough,  F.X.  Kneizys,  and  L.  S.  Rothman, 
FASCODE:  Fast  Atmospheric  Signature  Code  ( Spectral  Transmittance  and  Radiance ), 
AFGL-TR-78-008 1 ,  Air  Force  Geophysics  Laboratory,  Hanscom  AFB,  Massachusetts 
1978.  ADA057506. 

Stamnes,  K.,  S.C.  Tsay,  W.J.  Wiscombe,  and  K.  Jayaweerra,  Numerically  Stable  Algorithm 
for  Discrete-Ordinate-Method  Radiative  Transfer  in  Multiple  Scattering  and  Emitting 
Layered  Media,  Appl.  Opt.,  27,  2502-2509,  1988. 

Toon,  G.,  Jet  Propulsion  Laboratory,  Private  Communication,  November  -  December  1995. 

Varanasi,  P.,  and  V.  Nemtchinov,  Thermal  Infrared  Absorption  Coefficients  of  CFC-12  at 
Atmospheric  Conditions,  J.  Quant.  Spectrosc.  Radiat.  Transfer,  51,  679-687,  1994. 

Wang,  J.,  Phillips  Laboratory  Directorate  of  Geophysics,  Private  Communication,  1994. 


61 


62 


Appendix  A:  Spectral  Grid  Algorithm 


1°  FASCODE,  the  spectral  grid  is  optimized  individually  for  each  layer.  Layer  spectral 
quantities  (optical  depth,  transmittance  or  radiance)  are  combined  by  interpolating  the  coarser 
grid  onto  the  finer  grid.  The  layer  dv  is  set  by  default  to  one  fourth  the  average  Voigt  halfwidth 
for  that  layer.  This  value  represents  a  compromise  between  adequately  sampling  the  spectrum 
and  minimizing  the  number  of  spectral  points.  For  example,  for  the  US  Standard  Atmosphere, 
the  dv  varies  between  about  0.03  cm'1  at  the  ground  and  about  0.0003  cm’1  at  30  km. 

At  the  time  FASCODE  was  first  designed,  computer  memory  and  disk  space  were  limited 
and  minimizing  the  number  of  spectral  points  was  critical.  However,  the  cost  of  keeping  a 
separate  dv  for  each  layer  is  the  need  for  repeated  spectral  interpolation  when  combining  layers. 
Using  the  same  dv  for  all  layers  would  eliminate  the  need  for  spectral  interpolation,  a  viable 
alternative  since  the  increased  requirements  for  memory  and  disk  space  are  not  constraining 
factors  on  current  computer  systems. 

This  system  of  using  the  same  spectral  grid  spacing  for  all  layers  was  implemented  in  the 
model  XFWD,  a  modified  version  of  FASCODE  developed  internally  at  AER  (Miller  et  al., 
1999).  Both  models  calculate  the  optical  depth  at  a  spectral  resolution  appropriate  for  the  layer. 
XFWD  then  interpolates  these  values  to  the  final  resolution  before  computing  the  radiative 
transfer.  In  contrast,  FASCODE  computes  the  radiative  transfer  by  interpolating  the  results  from 
the  previous  layer  onto  the  grid  of  the  current  layer  and  merging  the  results.  The  result  of  this 
merge  is  written  to  a  "scratch"  file.  This  procedure  is  repeated  from  the  lowest  altitude  layer 
(largest  dv)  to  the  highest  (smallest  dv).  In  XFWD  all  of  the  optical  depth  calculations  are 
completed  before  the  radiative  transfer  is  done,  and  all  calculations  are  carried  out  in  memory. 
The  difference  in  doing  the  radiative  transfer  on  a  fixed  wavenumber  grid  can  be  seen  by 
comparing  the  times  to  compute  the  optical  depth  and  the  radiative  transfer. 

Timing  comparisons  were  conducted  between  FASCODE  and  XFWD  in  order  to  see  how 
the  proposed  program  structure  and  the  output  of  optical  depths  at  the  same  spectral  resolution 
for  each  layer  would  change  the  CPU  time.  The  timing  tests  were  conducted  for  calculations 
over  the  spectral  region  1010  -  1130  cm  \  To  eliminate  subtle  differences  in  the  methods  of 
computing  molecular  amounts,  the  same  atmosphere  was  input  to  both  models  (11  layers  for  the 
limb-viewing  case,  15  layers  for  the  up-  and  down-looking  cases,  all  with  U.S.  Standard 
Atmosphere  conditions)  and  the  same  line  file  was  used  (first  12  species  from  the  HITRAN 
database).  The  continuum  was  not  implemented,  but  in  both  cases  it  is  only  a  minor  contributor 
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to  the  total  computation  time.  The  cases  tested  are  representative  of  the  type  of  calculations  for 
which  FASCODE  is  used.  The  fact  that  not  all  of  the  features  were  implemented  should  not 
impact  the  timing  results;  while  only  a  single  comparison  is  discussed  for  each  case,  the  results 
agree  with  numerous  other  runs  of  the  two  models.  Table  A-  1  provides  the  total  CPU  time 
required  for  both  models  when  each  of  the  test  cases  was  run.  The  models  were  both  found  to 
use  the  same  percentage  of  CPU  time  relative  to  actual  time,  indicating  that  system  usage  did  not 
contribute  to  timing  differences.  Both  models  were  run  on  a  "SPARCcenter  1000"  computer. 
These  times  are  shown  for  each  of  the  test  case  conditions  in  Table  A-  2,  Table  A-  3,  and  Table 
A- 4. 


Table  A- 1  Model  Comparisons  for  Total  CPU  Run  Time 


CASE 

FASCOD3P  (seconds) 

XFWD(seconds) 

12  km  tangent  height 

275.53 

188.33 

0  km  - 100  km  (up-looking) 

309.22 

186.12 

100  km  -  0  km  (down-looking) 

317.03 

215.21 

Table  A-  2  Casel:  Limb  Calculation  with  12km  Tangent  Height 


LAYER 

FASCOD3P 

XFWD 

Radiative  Transfer 

Radiative  Transfer 

1 

2.70 

0.98 

4.05 

3.61 

2 

3.76 

1.70 

4.40 

3.52 

3 

6.92 

3.38 

5.83 

3.49 

4 

10.28 

5.58 

8.70 

3.55 

5 

21.47 

10.57 

14.08 

3.54 

6 

21.21 

11.71 

17.05 

3.44 

7 

20.99 

11.62 

17.80 

3.49 

8 

21.08 

11.64 

18.43 

3.51 

9 

20.90 

11.54 

19.22 

3.58 

10 

24.98 

13.62 

19.48 

3.58 

11 

25.19 

13.72 

20.14 

3.85 

Total 

179.48 

96.06 

149.18 

39.16 

Percentage 

65% 

35% 

79% 

21% 
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Table  A-  3  Case  2:  Looking  Up:  Okm  -  100km 


LAYER 

FASCOD3P 

XFWD 

Radiative  Transfer 

Optical  Depth 

Radiative  Transfer 

1 

1.96 

0.26 

3.60 

1.99 

2 

1.75 

0.34 

2.34 

1.94 

3 

1.90 

0.43 

2.36 

1.86 

4 

2.19 

0.68 

2.48 

1.88 

5 

3.06 

1.24 

2.93 

1.97 

6 

5.23 

2.42 

4.33 

1.94 

7 

7.64 

3.74 

5.70 

1.94 

8 

11.01 

5.92 

8.74 

1.95 

9 

16.30 

8.44 

13.93 

1.85 

10 

25.35 

13.11 

16.87 

1.83 

11 

25.12 

13.83 

17.65 

1.83 

12 

25.25 

14.43 

18.22 

1.76 

13 

25.40 

13.92 

19.17 

1.79 

14 

24.87 

13.83 

19.59 

1.78 

15 

25.63 

14.02 

19.88 

1.98 

Total 

202.66 

106.61 

157.79 

28.29 

Percentage 

65% 

35% 

85% 

15% 
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Table  A-  4  Case  3:  Looking  Down:  100km  -  0km 


LAYER 

FASCOD3P 

XFWD 

Radiative  Transfer 

Radiative  Transfer 

1 

2.15 

0.28 

4.15 

2.16 

2 

1.87 

0.35 

4.06 

2.12 

3 

2.08 

0.45 

4.21 

2.00 

4 

2.27 

0.67 

4.25 

2.06 

5 

3.22 

1.26 

4.74 

2.00 

6 

5.47 

2.60 

6.21 

1.99 

7 

7.74 

3.79 

7.69 

1.93 

8 

11.43 

5.78 

10.86 

1.89 

9 

16.79 

8.69 

15.97 

1.81 

10 

25.91 

13.31 

19.04 

1.78 

11 

25.51 

15.30 

20.08 

1.70 

12 

25.64 

14.16 

18.52 

1.65 

13 

25.38 

14.26 

21.34 

1.62 

14 

26.49 

14.25 

21.90 

1.58 

15 

25.87 

14.07 

24.12 

1.81 

Total 

207.82 

109.22 

187.14 

28.10 

66% 

34% 

87% 

13% 

The  timing  for  the  optical  depth  calculations  is  similar  in  FASCODE  and  XFWD  because 
they  are  calculated  at  the  same  monochromatic  spectral  resolution  (DV)  in  each  program.  Minor 
differences  in  the  timing  occur  because  the  optical  depths  are  interpolated  to  the  final  dv  in 
XFWD  while  they  are  output  to  disk  in  FASCODE.  However,  what  is  important  to  note  is  that 
for  a  particular  calculation  the  radiative  transfer  time  is  a  constant  for  XFWD  while  it  increases 
in  FASCODE  as  the  DV  decreases.  The  FASCODE  timing  for  radiative  transfer  listed  in  the 
above  tables  includes  I/O  time  while  in  the  case  of  XFWD  there  is  no  I/O.  Looking  at  Case  1,  in 
Table  A-  2,  we  can  assume  that  since  the  dv’s  for  the  two  models  are  the  same  for  Layer  1 1,  the 
time  to  compute  the  radiative  transfer  should  also  be  the  same.  This  implies  that  9.87  seconds  in 
FASCODE  (out  of  a  total  of  13.72  seconds)  can  be  attributed  to  the  I/Os.  Thus  it  is  clear  that  a 
significant  fraction  of  the  radiative  transfer  computation  time  in  FASCODE  consists  of  I/O, 
especially  at  lower  pressures. 

From  these  timing  results  we  may  make  the  following  conclusions: 
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(1)  The  I/Os  are  a  limiting  factor  and  make  up  more  than  80%  of  the  total  FASCODE 
radiative  transfer  time  in  up-  and  down-looking  cases. 

(2)  I/Os  are  suppressed  in  XFWD  by  making  optimal  use  of  the  available  core  memory. 
However,  the  radiative  transfer  time  in  XFWD  is  affected  by  the  fact  that  the  calculations  are 
done  at  high  spectral  resolution.  We  estimate  that  this  feature  accounts  for  an  increase  of 
about  100%  in  the  radiative  transfer  time  over  that  of  FASCODE  (not  including  FASCODE 
I/O).  Nonetheless  radiative  transfer  represents  only  a  small  fraction  (13-21%)  of  the  total 
computation  time,  which  indicates  that  the  impact  of  this  structure  on  the  total  time  is 
minimal  (7-10%). 

One  should  also  note  that  the  dv  is  pre-set  by  the  user  in  XFWD.  In  the  case  where  the 
atmosphere  is  truncated  at,  for  example,  Layer  2,  the  optical  depth  would  be  interpolated  to  a  dv 
appropriate  for  Layer  2  and  not  the  high  dv  of  Layer  1 5  as  in  the  above  examples.  In  general,  dv 
should  never  be  set  smaller  than  the  smallest  dv  in  the  layering  scheme. 

Comparisons  of  FASCODE  and  LBLRTM  timing  using  a  "benchmark"  input  developed  for 
LBLRTM  were  conducted  on  a  (non-vectorized)  "SPARCcenter  1000"  computer.  Tests  were 
also  done  to  compare  timing  differences  resulting  from  changes  to  the  Voigt  lineshape  algorithm. 
These  tests  indicated  that  changes  in  the  Voigt  algorithm  resulted  in  an  increase  of  about  10%  in 
the  LBLRTM  computation  time.  This  increase  was  subsequently  reduced  to  about  7%  through 
the  efficient  re-coding  of  some  sections  of  the  algorithm.  Improvements  in  the  structure  of  the 
LBLRTM  read  and  write  routines  ("BUFIN"  and  "BUFOUT")  should  also  favorably  impact  the 
timing  on  non-vectorized  machines.  In  general,  the  total  CPU  time  for  LBLRTM  (after  the  Voigt 
changes)  is  about  a  factor  of  two  less  than  that  for  an  identical  FASCODE  calculation.  These 
improvements  in  the  code  to  reduce  the  computation  time  are  included  in  FASE. 
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